//+------------------------------------------------------------------+
//|                                                  integration.mqh |
//|            Copyright 2003-2022 Sergey Bochkanov (ALGLIB project) |
//|                             Copyright 2012-2023, MetaQuotes Ltd. |
//|                                             https://www.mql5.com |
//+------------------------------------------------------------------+
//| Implementation of ALGLIB library in MetaQuotes Language 5        |
//|                                                                  |
//| The features of the library include:                             |
//| - Linear algebra (direct algorithms, EVD, SVD)                   |
//| - Solving systems of linear and non-linear equations             |
//| - Interpolation                                                  |
//| - Optimization                                                   |
//| - FFT (Fast Fourier Transform)                                   |
//| - Numerical integration                                          |
//| - Linear and nonlinear least-squares fitting                     |
//| - Ordinary differential equations                                |
//| - Computation of special functions                               |
//| - Descriptive statistics and hypothesis testing                  |
//| - Data analysis - classification, regression                     |
//| - Implementing linear algebra algorithms, interpolation, etc.    |
//|   in high-precision arithmetic (using MPFR)                      |
//|                                                                  |
//| This file is free software; you can redistribute it and/or       |
//| modify it under the terms of the GNU General Public License as   |
//| published by the Free Software Foundation (www.fsf.org); either  |
//| version 2 of the License, or (at your option) any later version. |
//|                                                                  |
//| This program is distributed in the hope that it will be useful,  |
//| but WITHOUT ANY WARRANTY; without even the implied warranty of   |
//| MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the     |
//| GNU General Public License for more details.                     |
//+------------------------------------------------------------------+
#include "ap.h"
#include "alglibinternal.h"
#include "linalg.h"
#include "specialfunctions.h"
//+------------------------------------------------------------------+
//| Gauss quadrature formula                                         |
//+------------------------------------------------------------------+
class CGaussQ
  {
public:
   static void       GQGenerateRec(double &alpha[],double &beta[],const double mu0,const int n,int &info,double &x[],double &w[]);
   static void       GQGenerateGaussLobattoRec(double &calpha[],double &cbeta[],const double mu0,const double a,const double b,int n,int &info,double &x[],double &w[]);
   static void       GQGenerateGaussRadauRec(double &calpha[],double &cbeta[],const double mu0,const double a,int n,int &info,double &x[],double &w[]);
   static void       GQGenerateGaussLegendre(const int n,int &info,double &x[],double &w[]);
   static void       GQGenerateGaussJacobi(const int n,const double alpha,const double beta,int &info,double &x[],double &w[]);
   static void       GQGenerateGaussLaguerre(const int n,const double alpha,int &info,double &x[],double &w[]);
   static void       GQGenerateGaussHermite(const int n,int &info,double &x[],double &w[]);
  };
//+------------------------------------------------------------------+
//| Computation of nodes and weights for a Gauss quadrature formula  |
//| The algorithm generates the N-point Gauss quadrature formula     |
//| with weight function given by coefficients alpha and beta of a   |
//| recurrence relation which generates a system of orthogonal       |
//| polynomials:                                                     |
//| P-1(x)   =  0                                                    |
//| P0(x)    =  1                                                    |
//| Pn+1(x)  =  (x-alpha(n))*Pn(x)  -  beta(n)*Pn-1(x)               |
//| and zeroth moment Mu0                                            |
//| Mu0 = integral(W(x)dx,a,b)                                       |
//| INPUT PARAMETERS:                                                |
//|     Alpha   -   array[0..N-1], alpha coefficients                |
//|     Beta    -   array[0..N-1], beta coefficients                 |
//|                 Zero-indexed element is not used and may be      |
//|                 arbitrary. Beta[I]>0.                            |
//|     Mu0     -   zeroth moment of the weight function.            |
//|     N       -   number of nodes of the quadrature formula, N>=1  |
//| OUTPUT PARAMETERS:                                               |
//|     Info    -   error code:                                      |
//|                 * -3    internal eigenproblem solver hasn't      |
//|                         converged                                |
//|                 * -2    Beta[i]<=0                               |
//|                 * -1    incorrect N was passed                   |
//|                 *  1    OK                                       |
//|     X       -   array[0..N-1] - array of quadrature nodes,       |
//|                 in ascending order.                              |
//|     W       -   array[0..N-1] - array of quadrature weights.     |
//+------------------------------------------------------------------+
void CGaussQ::GQGenerateRec(double &alpha[],double &beta[],
                            const double mu0,const int n,
                            int &info,double &x[],double &w[])
  {
//--- create a variable
   int i=0;
//--- create arrays
   double d[];
   double e[];
//--- create matrix
   CMatrixDouble z;
//--- initialization
   info=0;
//--- check
   if(n<1)
     {
      info=-1;
      return;
     }
   info=1;
//--- Initialize
   ArrayResize(d,n);
   ArrayResize(e,n);
//--- calculation
   for(i=1; i<=n-1; i++)
     {
      d[i-1]=alpha[i-1];
      //--- check
      if(beta[i]<=0.0)
        {
         info=-2;
         return;
        }
      e[i-1]=MathSqrt(beta[i]);
     }
//--- change value
   d[n-1]=alpha[n-1];
//--- EVD
   if(!CEigenVDetect::SMatrixTdEVD(d,e,n,3,z))
     {
      info=-3;
      return;
     }
//--- Generate
   ArrayResize(x,n);
   ArrayResize(w,n);
//--- calculation
   for(i=1; i<=n; i++)
     {
      x[i-1]=d[i-1];
      w[i-1]=mu0*CMath::Sqr(z[0][i-1]);
     }
  }
//+------------------------------------------------------------------+
//| Computation of nodes and weights for a Gauss-Lobatto quadrature  |
//| formula                                                          |
//| The algorithm generates the N-point Gauss-Lobatto quadrature     |
//| formula with weight function given by coefficients alpha and beta|
//| of a recurrence which generates a system of orthogonal           |
//| polynomials.                                                     |
//| P-1(x)   =  0                                                    |
//| P0(x)    =  1                                                    |
//| Pn+1(x)  =  (x-alpha(n))*Pn(x)  -  beta(n)*Pn-1(x)               |
//| and zeroth moment Mu0                                            |
//| Mu0 = integral(W(x)dx,a,b)                                       |
//| INPUT PARAMETERS:                                                |
//|     Alpha   -   array[0..N-2], alpha coefficients                |
//|     Beta    -   array[0..N-2], beta coefficients.                |
//|                 Zero-indexed element is not used, may be         |
//|                 arbitrary. Beta[I]>0                             |
//|     Mu0     -   zeroth moment of the weighting function.         |
//|     A       -   left boundary of the integration interval.       |
//|     B       -   right boundary of the integration interval.      |
//|     N       -   number of nodes of the quadrature formula, N>=3  |
//|                 (including the left and right boundary nodes).   |
//| OUTPUT PARAMETERS:                                               |
//|     Info    -   error code:                                      |
//|                 * -3    internal eigenproblem solver hasn't      |
//|                         converged                                |
//|                 * -2    Beta[i]<=0                               |
//|                 * -1    incorrect N was passed                   |
//|                 *  1    OK                                       |
//|     X       -   array[0..N-1] - array of quadrature nodes,       |
//|                 in ascending order.                              |
//|     W       -   array[0..N-1] - array of quadrature weights.     |
//+------------------------------------------------------------------+
void CGaussQ::GQGenerateGaussLobattoRec(double &calpha[],double &cbeta[],
                                        const double mu0,const double a,
                                        const double b,int n,int &info,
                                        double &x[],double &w[])
  {
//--- create variables
   int    i=0;
   double pim1a=0;
   double pia=0;
   double pim1b=0;
   double pib=0;
   double t=0;
   double a11=0;
   double a12=0;
   double a21=0;
   double a22=0;
   double b1=0;
   double b2=0;
   double alph=0;
   double bet=0;
//--- create arrays
   double d[];
   double e[];
//--- create matrix
   CMatrixDouble z;
//--- creating copies of arrays
   double alpha[];
   double beta[];
   ArrayCopy(alpha,calpha);
   ArrayCopy(beta,cbeta);
//--- initialization
   info=0;
//--- check
   if(n<=2)
     {
      info=-1;
      return;
     }
   info=1;
//--- Initialize,D[1:N+1],E[1:N]
   n=n-2;
//--- allocation
   ArrayResize(d,n+2);
   ArrayResize(e,n+1);
//--- copy
   for(i=1; i<=n+1; i++)
      d[i-1]=alpha[i-1];
   for(i=1; i<=n; i++)
     {
      //--- check
      if(beta[i]<=0.0)
        {
         info=-2;
         return;
        }
      e[i-1]=MathSqrt(beta[i]);
     }
//--- Caclulate Pn(a),Pn+1(a),Pn(b),Pn+1(b)
   beta[0]=0;
   pim1a=0;
   pia=1;
   pim1b=0;
   pib=1;
   for(i=1; i<=n+1; i++)
     {
      //--- Pi(a)
      t=(a-alpha[i-1])*pia-beta[i-1]*pim1a;
      pim1a=pia;
      pia=t;
      //--- Pi(b)
      t=(b-alpha[i-1])*pib-beta[i-1]*pim1b;
      pim1b=pib;
      pib=t;
     }
//--- Calculate alpha'(n+1),beta'(n+1)
   a11=pia;
   a12=pim1a;
   a21=pib;
   a22=pim1b;
   b1=a*pia;
   b2=b*pib;
//--- check
   if(MathAbs(a11)>MathAbs(a21))
     {
      //--- change values
      a22=a22-a12*a21/a11;
      b2=b2-b1*a21/a11;
      bet=b2/a22;
      alph=(b1-bet*a12)/a11;
     }
   else
     {
      //--- change values
      a12=a12-a22*a11/a21;
      b1=b1-b2*a11/a21;
      bet=b1/a12;
      alph=(b2-bet*a22)/a21;
     }
//--- check
   if(bet<0.0)
     {
      info=-3;
      return;
     }
//--- change values
   d[n+1]=alph;
   e[n]=MathSqrt(bet);
//--- EVD
   if(!CEigenVDetect::SMatrixTdEVD(d,e,n+2,3,z))
     {
      info=-3;
      return;
     }
//--- Generate
   ArrayResize(x,n+2);
   ArrayResize(w,n+2);
//--- get result
   for(i=1; i<=n+2; i++)
     {
      x[i-1]=d[i-1];
      w[i-1]=mu0*CMath::Sqr(z[0][i-1]);
     }
  }
//+------------------------------------------------------------------+
//| Computation of nodes and weights for a Gauss-Radau quadrature    |
//| formula                                                          |
//| The algorithm generates the N-point Gauss-Radau  quadrature      |
//| formula with weight function given by the coefficients alpha and |
//| beta of a recurrence which generates a system of orthogonal      |
//| polynomials.                                                     |
//| P-1(x)   =  0                                                    |
//| P0(x)    =  1                                                    |
//| Pn+1(x)  =  (x-alpha(n))*Pn(x)  -  beta(n)*Pn-1(x)               |
//| and zeroth moment Mu0                                            |
//| Mu0 = integral(W(x)dx,a,b)                                       |
//| INPUT PARAMETERS:                                                |
//|     Alpha   -   array[0..N-2], alpha coefficients.               |
//|     Beta    -   array[0..N-1], beta coefficients                 |
//|                 Zero-indexed element is not used.                |
//|                 Beta[I]>0                                        |
//|     Mu0     -   zeroth moment of the weighting function.         |
//|     A       -   left boundary of the integration interval.       |
//|     N       -   number of nodes of the quadrature formula, N>=2  |
//|                 (including the left boundary node).              |
//| OUTPUT PARAMETERS:                                               |
//|     Info    -   error code:                                      |
//|                 * -3    internal eigenproblem solver hasn't      |
//|                         converged                                |
//|                 * -2    Beta[i]<=0                               |
//|                 * -1    incorrect N was passed                   |
//|                 *  1    OK                                       |
//|     X       -   array[0..N-1] - array of quadrature nodes,       |
//|                 in ascending order.                              |
//|     W       -   array[0..N-1] - array of quadrature weights.     |
//+------------------------------------------------------------------+
void CGaussQ::GQGenerateGaussRadauRec(double &calpha[],double &cbeta[],
                                      const double mu0,const double a,
                                      int n,int &info,double &x[],
                                      double &w[])
  {
//--- create variables
   int i=0;
   double polim1=0;
   double poli=0;
   double t=0;
//--- create arrays
   double d[];
   double e[];
//--- create matrix
   CMatrixDouble z;
//--- creating copies of arrays
   double alpha[];
   double beta[];
   ArrayCopy(alpha,calpha);
   ArrayCopy(beta,cbeta);
//--- initialization
   info=0;
//--- check
   if(n<2)
     {
      info=-1;
      return;
     }
   info=1;
//--- Initialize,D[1:N],E[1:N]
   n=n-1;
//--- allocation
   ArrayResize(d,n+1);
   ArrayResize(e,n);
   for(i=1; i<=n; i++)
     {
      d[i-1]=alpha[i-1];
      //--- check
      if(beta[i]<=0.0)
        {
         info=-2;
         return;
        }
      e[i-1]=MathSqrt(beta[i]);
     }
//--- Caclulate Pn(a),Pn-1(a),and D[N+1]
   beta[0]=0;
   polim1=0;
   poli=1;
//--- calculation
   for(i=1; i<=n; i++)
     {
      t=(a-alpha[i-1])*poli-beta[i-1]*polim1;
      polim1=poli;
      poli=t;
     }
   d[n]=a-beta[n]*polim1/poli;
//--- EVD
   if(!CEigenVDetect::SMatrixTdEVD(d,e,n+1,3,z))
     {
      info=-3;
      return;
     }
//--- Generate
   ArrayResize(x,n+1);
   ArrayResize(w,n+1);
//--- get result
   for(i=1; i<=n+1; i++)
     {
      x[i-1]=d[i-1];
      w[i-1]=mu0*CMath::Sqr(z[0][i-1]);
     }
  }
//+------------------------------------------------------------------+
//| Returns nodes/weights for Gauss-Legendre quadrature on [-1,1]    |
//| with N nodes.                                                    |
//| INPUT PARAMETERS:                                                |
//|     N           -   number of nodes, >=1                         |
//| OUTPUT PARAMETERS:                                               |
//|     Info        -   error code:                                  |
//|                     * -4    an error was detected when           |
//|                             calculating weights/nodes. N is too  |
//|                             large to obtain weights/nodes with   |
//|                             high enough accuracy. Try to use     |
//|                             multiple precision version.          |
//|                     * -3    internal eigenproblem solver hasn't  |
//|                             converged                            |
//|                     * -1    incorrect N was passed               |
//|                     * +1    OK                                   |
//|     X           -   array[0..N-1] - array of quadrature nodes,   |
//|                     in ascending order.                          |
//|     W           -   array[0..N-1] - array of quadrature weights. |
//+------------------------------------------------------------------+
void CGaussQ::GQGenerateGaussLegendre(const int n,int &info,
                                      double &x[],double &w[])
  {
//--- create a variable
   int i=0;
//--- create arrays
   double alpha[];
   double beta[];
//--- initialization
   info=0;
//--- check
   if(n<1)
     {
      info=-1;
      return;
     }
//--- allocation
   ArrayResize(alpha,n);
   ArrayResize(beta,n);
//--- initialization
   for(i=0; i<=n-1; i++)
      alpha[i]=0;
   beta[0]=2;
   for(i=1; i<=n-1; i++)
      beta[i]=1/(4-1/CMath::Sqr(i));
//--- function call
   GQGenerateRec(alpha,beta,beta[0],n,info,x,w);
//--- test basic properties to detect errors
   if(info>0)
     {
      //--- check
      if(x[0]<-1.0 || x[n-1]>1.0)
         info=-4;
      for(i=0; i<=n-2; i++)
        {
         //--- check
         if(x[i]>=x[i+1])
            info=-4;
        }
     }
  }
//+------------------------------------------------------------------+
//| Returns  nodes/weights  for  Gauss-Jacobi quadrature on [-1,1]   |
//| with weight function W(x)=Power(1-x,Alpha)*Power(1+x,Beta).      |
//| INPUT PARAMETERS:                                                |
//|     N           -   number of nodes, >=1                         |
//|     Alpha       -   power-law coefficient, Alpha>-1              |
//|     Beta        -   power-law coefficient, Beta>-1               |
//| OUTPUT PARAMETERS:                                               |
//|     Info        -   error code:                                  |
//|                     * -4    an error was detected when           |
//|                             calculating weights/nodes. Alpha or  |
//|                             Beta are too close to -1 to obtain   |
//|                             weights/nodes with high enough       |
//|                             accuracy, or, may be, N is too large.|
//|                             Try to use multiple precision version|
//|                     * -3    internal eigenproblem solver hasn't  |
//|                             converged                            |
//|                     * -1    incorrect N/Alpha/Beta was passed    |
//|                     * +1    OK                                   |
//|     X           -   array[0..N-1] - array of quadrature nodes,   |
//|                     in ascending order.                          |
//|     W           -   array[0..N-1] - array of quadrature weights. |
//+------------------------------------------------------------------+
void CGaussQ::GQGenerateGaussJacobi(const int n,const double alpha,
                                    const double beta,int &info,
                                    double &x[],double &w[])
  {
//--- create variables
   double alpha2=0;
   double beta2=0;
   double apb=0;
   double t=0;
   int    i=0;
   double s=0;
//--- create arrays
   double a[];
   double b[];
//--- initialization
   info=0;
//--- check
   if((n<1 || alpha<=-1.0) || beta<=-1.0)
     {
      info=-1;
      return;
     }
//--- allocation
   ArrayResize(a,n);
   ArrayResize(b,n);
//--- initialization
   apb=alpha+beta;
   a[0]=(beta-alpha)/(apb+2);
   t=(apb+1)*MathLog(2)+CGammaFunc::LnGamma(alpha+1,s)+CGammaFunc::LnGamma(beta+1,s)-CGammaFunc::LnGamma(apb+2,s);
//--- check
   if(t>MathLog(CMath::m_maxrealnumber))
     {
      info=-4;
      return;
     }
   b[0]=MathExp(t);
//--- check
   if(n>1)
     {
      //--- change values
      alpha2=CMath::Sqr(alpha);
      beta2=CMath::Sqr(beta);
      a[1]=(beta2-alpha2)/((apb+2)*(apb+4));
      b[1]=4*(alpha+1)*(beta+1)/((apb+3)*CMath::Sqr(apb+2));
      //--- calculation
      for(i=2; i<=n-1; i++)
        {
         a[i]=0.25*(beta2-alpha2)/(i*i*(1+0.5*apb/i)*(1+0.5*(apb+2)/i));
         b[i]=0.25*(1+alpha/i)*(1+beta/i)*(1+apb/i)/((1+0.5*(apb+1)/i)*(1+0.5*(apb-1)/i)*CMath::Sqr(1+0.5*apb/i));
        }
     }
//--- function call
   GQGenerateRec(a,b,b[0],n,info,x,w);
//--- test basic properties to detect errors
   if(info>0)
     {
      //--- check
      if(x[0]<-1.0 || x[n-1]>1.0)
         info=-4;
      for(i=0; i<=n-2; i++)
        {
         //--- check
         if(x[i]>=x[i+1])
            info=-4;
        }
     }
  }
//+------------------------------------------------------------------+
//| Returns nodes/weights for Gauss-Laguerre quadrature on [0,+inf)  |
//| with weight function W(x)=Power(x,Alpha)*Exp(-x)                 |
//| INPUT PARAMETERS:                                                |
//|     N           -   number of nodes, >=1                         |
//|     Alpha       -   power-law coefficient, Alpha>-1              |
//| OUTPUT PARAMETERS:                                               |
//|     Info        -   error code:                                  |
//|                     * -4    an error was detected when           |
//|                             calculating weights/nodes. Alpha is  |
//|                             too close to -1 to obtain            |
//|                             weights/nodes with high enough       |
//|                             accuracy or, may  be, N is too large.|
//|                             Try  to  use multiple precision      |
//|                             version.                             |
//|                     * -3    internal eigenproblem solver hasn't  |
//|                             converged                            |
//|                     * -1    incorrect N/Alpha was passed         |
//|                     * +1    OK                                   |
//|     X           -   array[0..N-1] - array of quadrature nodes,   |
//|                     in ascending order.                          |
//|     W           -   array[0..N-1] - array of quadrature weights. |
//+------------------------------------------------------------------+
void CGaussQ::GQGenerateGaussLaguerre(const int n,const double alpha,
                                      int &info,double &x[],
                                      double &w[])
  {
//--- create arrays
   double a[];
   double b[];
//--- create variables
   double t=0;
   int    i=0;
   double s=0;
//--- initialization
   info=0;
//--- check
   if(n<1 || alpha<=-1.0)
     {
      info=-1;
      return;
     }
//--- allocation
   ArrayResize(a,n);
   ArrayResize(b,n);
//--- initialization
   a[0]=alpha+1;
   t=CGammaFunc::LnGamma(alpha+1,s);
//--- check
   if(t>=MathLog(CMath::m_maxrealnumber))
     {
      info=-4;
      return;
     }
   b[0]=MathExp(t);
//--- check
   if(n>1)
     {
      //--- calculation
      for(i=1; i<=n-1; i++)
        {
         a[i]=2*i+alpha+1;
         b[i]=i*(i+alpha);
        }
     }
//--- function call
   GQGenerateRec(a,b,b[0],n,info,x,w);
//--- test basic properties to detect errors
   if(info>0)
     {
      //--- check
      if(x[0]<0.0)
         info=-4;
      for(i=0; i<=n-2; i++)
        {
         //--- check
         if(x[i]>=x[i+1])
            info=-4;
        }
     }
  }
//+------------------------------------------------------------------+
//| Returns  nodes/weights  for  Gauss-Hermite  quadrature on        |
//| (-inf,+inf) with weight function W(x)=Exp(-x*x)                  |
//| INPUT PARAMETERS:                                                |
//|     N           -   number of nodes, >=1                         |
//| OUTPUT PARAMETERS:                                               |
//|     Info        -   error code:                                  |
//|                     * -4    an error was detected when           |
//|                             calculating weights/nodes. May be, N |
//|                             is too large. Try to use multiple    |
//|                             precision version.                   |
//|                     * -3    internal eigenproblem solver hasn't  |
//|                             converged                            |
//|                     * -1    incorrect N/Alpha was passed         |
//|                     * +1    OK                                   |
//|     X           -   array[0..N-1] - array of quadrature nodes,   |
//|                     in ascending order.                          |
//|     W           -   array[0..N-1] - array of quadrature weights. |
//+------------------------------------------------------------------+
void CGaussQ::GQGenerateGaussHermite(const int n,int &info,
                                     double &x[],double &w[])
  {
//--- create a variable
   int i=0;
//--- create arrays
   double a[];
   double b[];
//--- initialization
   info=0;
//--- check
   if(n<1)
     {
      info=-1;
      return;
     }
//--- allocation
   ArrayResize(a,n);
   ArrayResize(b,n);
//--- initialization
   for(i=0; i<=n-1; i++)
      a[i]=0;
   b[0]=MathSqrt(4*MathArctan(1));
//--- check
   if(n>1)
     {
      for(i=1; i<=n-1; i++)
         b[i]=0.5*i;
     }
//--- function call
   GQGenerateRec(a,b,b[0],n,info,x,w);
//--- test basic properties to detect errors
   if(info>0)
     {
      for(i=0; i<=n-2; i++)
        {
         //--- check
         if(x[i]>=x[i+1])
            info=-4;
        }
     }
  }
//+------------------------------------------------------------------+
//| Gauss-Kronrod quadrature formula                                 |
//+------------------------------------------------------------------+
class CGaussKronrodQ
  {
public:
   static void       GKQGenerateRec(double &calpha[],double &cbeta[],const double mu0,int n,int &info,double &x[],double &wkronrod[],double &wgauss[]);
   static void       GKQGenerateGaussLegendre(const int n,int &info,double &x[],double &wkronrod[],double &wgauss[]);
   static void       GKQGenerateGaussJacobi(const int n,const double alpha,const double beta,int &info,double &x[],double &wkronrod[],double &wgauss[]);
   static void       GKQLegendreCalc(const int n,int &info,double &x[],double &wkronrod[],double &wgauss[]);
   static void       GKQLegendreTbl(const int n,double &x[],double &wkronrod[],double &wgauss[],double &eps);
  };
//+------------------------------------------------------------------+
//| Computation of nodes and weights of a Gauss-Kronrod quadrature   |
//| formula                                                          |
//| The algorithm generates the N-point Gauss-Kronrod quadrature     |
//| formula with weight function given by coefficients alpha and beta|
//| of a recurrence relation which generates a system of orthogonal  |
//| polynomials:                                                     |
//|     P-1(x)   =  0                                                |
//|     P0(x)    =  1                                                |
//|     Pn+1(x)  =  (x-alpha(n))*Pn(x)  -  beta(n)*Pn-1(x)           |
//| and zero moment Mu0                                              |
//|     Mu0 = integral(W(x)dx,a,b)                                   |
//| INPUT PARAMETERS:                                                |
//|     Alpha       -   alpha coefficients, array[0..floor(3*K/2)].  |
//|     Beta        -   beta coefficients,  array[0..ceil(3*K/2)].   |
//|                     Beta[0] is not used and may be arbitrary.    |
//|                     Beta[I]>0.                                   |
//|     Mu0         -   zeroth moment of the weight function.        |
//|     N           -   number of nodes of the Gauss-Kronrod         |
//|                     quadrature formula,                          |
//|                     N >= 3,                                      |
//|                     N =  2*K+1.                                  |
//| OUTPUT PARAMETERS:                                               |
//|     Info        -   error code:                                  |
//|                     * -5    no real and positive Gauss-Kronrod   |
//|                             formula can be created for such a    |
//|                             weight function with a given number  |
//|                             of nodes.                            |
//|                     * -4    N is too large, task may be ill      |
//|                             conditioned - x[i]=x[i+1] found.     |
//|                     * -3    internal eigenproblem solver hasn't  |
//|                             converged                            |
//|                     * -2    Beta[i]<=0                           |
//|                     * -1    incorrect N was passed               |
//|                     * +1    OK                                   |
//|     X           -   array[0..N-1] - array of quadrature nodes,   |
//|                     in ascending order.                          |
//|     WKronrod    -   array[0..N-1] - Kronrod weights              |
//|     WGauss      -   array[0..N-1] - Gauss weights (interleaved   |
//|                     with zeros corresponding to extended Kronrod |
//|                     nodes).                                      |
//+------------------------------------------------------------------+
void CGaussKronrodQ::GKQGenerateRec(double &calpha[],double &cbeta[],
                                    const double mu0,int n,int &info,
                                    double &x[],double &wkronrod[],
                                    double &wgauss[])
  {
//--- create variables
   int    i=0;
   int    j=0;
   int    wlen=0;
   int    woffs=0;
   double u=0;
   int    m=0;
   int    l=0;
   int    k=0;
   int    i_=0;
//--- create arrays
   double ta[];
   double t[];
   double s[];
   double xgtmp[];
   double wgtmp[];
//--- creating copies of arrays
   double alpha[];
   double beta[];
   ArrayCopy(alpha,calpha);
   ArrayCopy(beta,cbeta);
//--- initialization
   info=0;
//--- check
   if(n%2!=1 || n<3)
     {
      info=-1;
      return;
     }
   for(i=0; i<=(int)MathCeil((double)(3*(n/2))/2.0); i++)
     {
      //--- check
      if(beta[i]<=0.0)
        {
         info=-2;
         return;
        }
     }
   info=1;
//--- from external conventions about N/Beta/Mu0 to internal
   n=n/2;
   beta[0]=mu0;
//--- Calculate Gauss nodes/weights,save them for later processing
   CGaussQ::GQGenerateRec(alpha,beta,mu0,n,info,xgtmp,wgtmp);
//--- check
   if(info<0)
      return;
//--- Resize:
//--- * A from 0..floor(3*n/2) to 0..2*n
//--- * B from 0..ceil(3*n/2)  to 0..2*n
   ArrayResize(ta,(int)MathFloor((double)(3*n)/2.0)+1);
//--- copy
   for(i_=0; i_<=(int)MathFloor((double)(3*n)/2.0); i_++)
      ta[i_]=alpha[i_];
//--- allocation
   ArrayResize(alpha,2*n+1);
//--- copy
   for(i_=0; i_<=(int)MathFloor((double)(3*n)/2.0); i_++)
      alpha[i_]=ta[i_];
//--- initialization
   for(i=(int)MathFloor((double)(3*n)/2.0)+1; i<=2*n; i++)
      alpha[i]=0;
//--- allocation
   ArrayResize(ta,(int)MathCeil((double)(3*n)/2.0)+1);
//--- copy
   for(i_=0; i_<=(int)MathCeil((double)(3*n)/2.0); i_++)
      ta[i_]=beta[i_];
//--- allocation
   ArrayResize(beta,2*n+1);
//--- copy
   for(i_=0; i_<=(int)MathCeil((double)(3*n)/2.0); i_++)
      beta[i_]=ta[i_];
//--- initialization
   for(i=(int)MathCeil((double)(3*n)/2.0)+1; i<=2*n; i++)
      beta[i]=0;
//--- Initialize T,S
   wlen=2+n/2;
//--- allocation
   ArrayResize(t,wlen);
   ArrayResize(s,wlen);
   ArrayResize(ta,wlen);
   woffs=1;
   for(i=0; i<=wlen-1; i++)
     {
      t[i]=0;
      s[i]=0;
     }
//--- Algorithm from Dirk P. Laurie,"Calculation of Gauss-Kronrod quadrature rules",1997.
   t[woffs+0]=beta[n+1];
   for(m=0; m<=n-2; m++)
     {
      u=0;
      //--- calculation
      for(k=(m+1)/2; k>=0; k--)
        {
         l=m-k;
         u=u+(alpha[k+n+1]-alpha[l])*t[woffs+k]+beta[k+n+1]*s[woffs+k-1]-beta[l]*s[woffs+k];
         s[woffs+k]=u;
        }
      //--- copy
      for(i_=0; i_<=wlen-1; i_++)
         ta[i_]=t[i_];
      //--- copy
      for(i_=0; i_<=wlen-1; i_++)
         t[i_]=s[i_];
      //--- copy
      for(i_=0; i_<=wlen-1; i_++)
         s[i_]=ta[i_];
     }
//--- copy
   for(j=n/2; j>=0; j--)
      s[woffs+j]=s[woffs+j-1];
//--- cycle
   for(m=n-1; m<=2*n-3; m++)
     {
      u=0;
      //--- calculation
      for(k=m+1-n; k<=(m-1)/2; k++)
        {
         l=m-k;
         j=n-1-l;
         u=u-(alpha[k+n+1]-alpha[l])*t[woffs+j]-beta[k+n+1]*s[woffs+j]+beta[l]*s[woffs+j+1];
         s[woffs+j]=u;
        }
      //--- check
      if(m%2==0)
        {
         k=m/2;
         alpha[k+n+1]=alpha[k]+(s[woffs+j]-beta[k+n+1]*s[woffs+j+1])/t[woffs+j+1];
        }
      else
        {
         k=(m+1)/2;
         beta[k+n+1]=s[woffs+j]/s[woffs+j+1];
        }
      //--- copy
      for(i_=0; i_<=wlen-1; i_++)
         ta[i_]=t[i_];
      //--- copy
      for(i_=0; i_<=wlen-1; i_++)
         t[i_]=s[i_];
      //--- copy
      for(i_=0; i_<=wlen-1; i_++)
         s[i_]=ta[i_];
     }
//--- change value
   alpha[2*n]=alpha[n-1]-beta[2*n]*s[woffs+0]/t[woffs+0];
//--- calculation of Kronrod nodes and weights,unpacking of Gauss weights
   CGaussQ::GQGenerateRec(alpha,beta,mu0,2*n+1,info,x,wkronrod);
//--- check
   if(info==-2)
      info=-5;
//--- check
   if(info<0)
      return;
   for(i=0; i<=2*n-1; i++)
     {
      //--- check
      if(x[i]>=x[i+1])
         info=-4;
     }
//--- check
   if(info<0)
      return;
//--- allocation
   ArrayResize(wgauss,2*n+1);
//--- get result
   for(i=0; i<=2*n; i++)
      wgauss[i]=0;
   for(i=0; i<=n-1; i++)
      wgauss[2*i+1]=wgtmp[i];
  }
//+------------------------------------------------------------------+
//| Returns Gauss and Gauss-Kronrod nodes/weights for Gauss-Legendre |
//| quadrature with N points.                                        |
//| GKQLegendreCalc (calculation) or GKQLegendreTbl (precomputed     |
//| table) is used depending on machine precision and number of      |
//| nodes.                                                           |
//| INPUT PARAMETERS:                                                |
//|     N           -   number of Kronrod nodes, must be odd number, |
//|                     >=3.                                         |
//| OUTPUT PARAMETERS:                                               |
//|     Info        -   error code:                                  |
//|                     * -4    an error was detected when           |
//|                             calculating weights/nodes. N is too  |
//|                             obtain large to weights/nodes with   |
//|                             high enough accuracy. Try to use     |
//|                             multiple precision version.          |
//|                     * -3    internal eigenproblem solver hasn't  |
//|                             converged                            |
//|                     * -1    incorrect N was passed               |
//|                     * +1    OK                                   |
//|     X           -   array[0..N-1] - array of quadrature nodes,   |
//|                     ordered in ascending order.                  |
//|     WKronrod    -   array[0..N-1] - Kronrod weights              |
//|     WGauss      -   array[0..N-1] - Gauss weights (interleaved   |
//|                     with zeros corresponding to extended Kronrod |
//|                     nodes).                                      |
//+------------------------------------------------------------------+
void CGaussKronrodQ::GKQGenerateGaussLegendre(const int n,int &info,
                                              double &x[],double &wkronrod[],
                                              double &wgauss[])
  {
//--- create a variable
   double eps=0;
//--- initialization
   info=0;
//--- check
   if(CMath::m_machineepsilon>1.0E-32 && (((((n==15 || n==21) || n==31) || n==41) || n==51) || n==61))
     {
      info=1;
      GKQLegendreTbl(n,x,wkronrod,wgauss,eps);
     }
   else
      GKQLegendreCalc(n,info,x,wkronrod,wgauss);
  }
//+------------------------------------------------------------------+
//| Returns Gauss and Gauss-Kronrod nodes/weights for Gauss-Jacobi   |
//| quadrature on [-1,1] with weight function                        |
//|     W(x)=Power(1-x,Alpha)*Power(1+x,Beta).                       |
//| INPUT PARAMETERS:                                                |
//|     N           -   number of Kronrod nodes, must be odd number, |
//|                     >=3.                                         |
//|     Alpha       -   power-law coefficient, Alpha>-1              |
//|     Beta        -   power-law coefficient, Beta>-1               |
//| OUTPUT PARAMETERS:                                               |
//|     Info        -   error code:                                  |
//|                     * -5    no real and positive Gauss-Kronrod   |
//|                             formula can be created for such a    |
//|                             weight function with a given number  |
//|                             of nodes.                            |
//|                     * -4    an  error was detected when          |
//|                             calculating weights/nodes. Alpha or  |
//|                             Beta are too close to -1 to obtain   |
//|                             weights/nodes with high enough       |
//|                             accuracy, or, may be, N is too large.|
//|                             Try to use multiple precision version|
//|                     * -3    internal eigenproblem solver hasn't  |
//|                             converged                            |
//|                     * -1    incorrect N was passed               |
//|                     * +1    OK                                   |
//|                     * +2    OK, but quadrature rule have exterior|
//|                             nodes, x[0]<-1 or x[n-1]>+1          |
//|     X           -   array[0..N-1] - array of quadrature nodes,   |
//|                     ordered in ascending order.                  |
//|     WKronrod    -   array[0..N-1] - Kronrod weights              |
//|     WGauss      -   array[0..N-1] - Gauss weights (interleaved   |
//|                     with zeros corresponding to extended Kronrod |
//|                     nodes).                                      |
//+------------------------------------------------------------------+
void CGaussKronrodQ::GKQGenerateGaussJacobi(const int n,const double alpha,
                                            const double beta,int &info,
                                            double &x[],double &wkronrod[],
                                            double &wgauss[])
  {
//--- create variables
   int    clen=0;
   double alpha2=0;
   double beta2=0;
   double apb=0;
   double t=0;
   int    i=0;
   double s=0;
//--- create arrays
   double a[];
   double b[];
//--- initialization
   info=0;
//--- check
   if(n%2!=1 || n<3)
     {
      info=-1;
      return;
     }
//--- check
   if(alpha<=-1.0 || beta<=-1.0)
     {
      info=-1;
      return;
     }
//--- change value
   clen=(int)MathCeil((double)(3*(n/2))/2.0)+1;
//--- allocation
   ArrayResize(a,clen);
   ArrayResize(b,clen);
//--- initialization
   for(i=0; i<=clen-1; i++)
      a[i]=0;
   apb=alpha+beta;
   a[0]=(beta-alpha)/(apb+2);
   t=(apb+1)*MathLog(2)+CGammaFunc::LnGamma(alpha+1,s)+CGammaFunc::LnGamma(beta+1,s)-CGammaFunc::LnGamma(apb+2,s);
//--- check
   if(t>MathLog(CMath::m_maxrealnumber))
     {
      info=-4;
      return;
     }
   b[0]=MathExp(t);
//--- check
   if(clen>1)
     {
      //--- change values
      alpha2=CMath::Sqr(alpha);
      beta2=CMath::Sqr(beta);
      a[1]=(beta2-alpha2)/((apb+2)*(apb+4));
      b[1]=4*(alpha+1)*(beta+1)/((apb+3)*CMath::Sqr(apb+2));
      //--- calculation
      for(i=2; i<=clen-1; i++)
        {
         a[i]=0.25*(beta2-alpha2)/(i*i*(1+0.5*apb/i)*(1+0.5*(apb+2)/i));
         b[i]=0.25*(1+alpha/i)*(1+beta/i)*(1+apb/i)/((1+0.5*(apb+1)/i)*(1+0.5*(apb-1)/i)*CMath::Sqr(1+0.5*apb/i));
        }
     }
//--- function call
   GKQGenerateRec(a,b,b[0],n,info,x,wkronrod,wgauss);
//--- test basic properties to detect errors
   if(info>0)
     {
      //--- check
      if(x[0]<-1.0 || x[n-1]>1.0)
         info=2;
      for(i=0; i<=n-2; i++)
        {
         //--- check
         if(x[i]>=x[i+1])
            info=-4;
        }
     }
  }
//+------------------------------------------------------------------+
//| Returns Gauss and Gauss-Kronrod nodes for quadrature with N      |
//| points.                                                          |
//| Reduction to tridiagonal eigenproblem is used.                   |
//| INPUT PARAMETERS:                                                |
//|     N           -   number of Kronrod nodes, must be odd number, |
//|                     >=3.                                         |
//| OUTPUT PARAMETERS:                                               |
//|     Info        -   error code:                                  |
//|                     * -4    an error was detected when           |
//|                             calculating weights/nodes. N is too  |
//|                             large to obtain weights/nodes with   |
//|                             high enough accuracy.                |
//|                             Try to use multiple precision        |
//|                             version.                             |
//|                     * -3    internal eigenproblem solver hasn't  |
//|                             converged                            |
//|                     * -1    incorrect N was passed               |
//|                     * +1    OK                                   |
//|     X           -   array[0..N-1] - array of quadrature nodes,   |
//|                     ordered in ascending order.                  |
//|     WKronrod    -   array[0..N-1] - Kronrod weights              |
//|     WGauss      -   array[0..N-1] - Gauss weights (interleaved   |
//|                     with zeros corresponding to extended Kronrod |
//|                     nodes).                                      |
//+------------------------------------------------------------------+
void CGaussKronrodQ::GKQLegendreCalc(const int n,int &info,double &x[],
                                     double &wkronrod[],double &wgauss[])
  {
//--- create variables
   int    alen=0;
   int    blen=0;
   double mu0=0;
   int    k=0;
   int    i=0;
//--- create arrays
   double alpha[];
   double beta[];
//--- initialization
   info=0;
//--- check
   if(n%2!=1 || n<3)
     {
      info=-1;
      return;
     }
//--- initialization
   mu0=2;
   alen=(int)MathFloor((double)(3*(n/2))/2.0)+1;
   blen=(int)MathCeil((double)(3*(n/2))/2.0)+1;
//--- allocation
   ArrayResize(alpha,alen);
   ArrayResize(beta,blen);
//--- initialization
   for(k=0; k<=alen-1; k++)
      alpha[k]=0;
   beta[0]=2;
   for(k=1; k<=blen-1; k++)
      beta[k]=1/(4-1/CMath::Sqr(k));
//--- function call
   GKQGenerateRec(alpha,beta,mu0,n,info,x,wkronrod,wgauss);
//--- test basic properties to detect errors
   if(info>0)
     {
      //--- check
      if(x[0]<-1.0 || x[n-1]>1.0)
         info=-4;
      for(i=0; i<=n-2; i++)
        {
         //--- check
         if(x[i]>=x[i+1])
            info=-4;
        }
     }
  }
//+------------------------------------------------------------------+
//| Returns Gauss and Gauss-Kronrod nodes for quadrature with N      |
//| points using pre-calculated table. Nodes/weights were computed   |
//| with accuracy up to 1.0E-32 (if MPFR version of ALGLIB is used). |
//| In standard double precision accuracy reduces to something about |
//| 2.0E-16 (depending  on your compiler's handling of long floating |
//| point constants).                                                |
//| INPUT PARAMETERS:                                                |
//|     N           -   number of Kronrod nodes.                     |
//|                     N can be 15, 21, 31, 41, 51, 61.             |
//| OUTPUT PARAMETERS:                                               |
//|     X           -   array[0..N-1] - array of quadrature nodes,   |
//|                     ordered in ascending order.                  |
//|     WKronrod    -   array[0..N-1] - Kronrod weights              |
//|     WGauss      -   array[0..N-1] - Gauss weights (interleaved   |
//|                     with zeros corresponding to extended Kronrod |
//|                     nodes).                                      |
//+------------------------------------------------------------------+
void CGaussKronrodQ::GKQLegendreTbl(const int n,double &x[],
                                    double &wkronrod[],double &wgauss[],
                                    double &eps)
  {
//--- create variables
   int    i=0;
   int    ng=0;
   double tmp=0;
//--- create arrays
   int p1[];
   int p2[];
//--- initialization
   eps=0;
//--- these initializers are not really necessary,
//--- but without them compiler complains about uninitialized locals
   ng=0;
//--- Process
   if(!CAp::Assert(n==15 || n==21 || n==31 || n==41 || n==51 || n==61,__FUNCTION__+": incorrect N!"))
      return;
//--- allocation
   ArrayResize(x,n);
   ArrayResize(wkronrod,n);
   ArrayResize(wgauss,n);
//--- initialization
   for(i=0; i<=n-1; i++)
     {
      x[i]=0;
      wkronrod[i]=0;
      wgauss[i]=0;
     }
   eps=MathMax(CMath::m_machineepsilon,1.0E-32);
//--- check
   if(n==15)
     {
      //--- change values
      ng=4;
      wgauss[0]=0.129484966168869693270611432679082;
      wgauss[1]=0.279705391489276667901467771423780;
      wgauss[2]=0.381830050505118944950369775488975;
      wgauss[3]=0.417959183673469387755102040816327;
      x[0]=0.991455371120812639206854697526329;
      x[1]=0.949107912342758524526189684047851;
      x[2]=0.864864423359769072789712788640926;
      x[3]=0.741531185599394439863864773280788;
      x[4]=0.586087235467691130294144838258730;
      x[5]=0.405845151377397166906606412076961;
      x[6]=0.207784955007898467600689403773245;
      x[7]=0.000000000000000000000000000000000;
      wkronrod[0]=0.022935322010529224963732008058970;
      wkronrod[1]=0.063092092629978553290700663189204;
      wkronrod[2]=0.104790010322250183839876322541518;
      wkronrod[3]=0.140653259715525918745189590510238;
      wkronrod[4]=0.169004726639267902826583426598550;
      wkronrod[5]=0.190350578064785409913256402421014;
      wkronrod[6]=0.204432940075298892414161999234649;
      wkronrod[7]=0.209482141084727828012999174891714;
     }
//--- check
   if(n==21)
     {
      //--- change values
      ng=5;
      wgauss[0]=0.066671344308688137593568809893332;
      wgauss[1]=0.149451349150580593145776339657697;
      wgauss[2]=0.219086362515982043995534934228163;
      wgauss[3]=0.269266719309996355091226921569469;
      wgauss[4]=0.295524224714752870173892994651338;
      x[0]=0.995657163025808080735527280689003;
      x[1]=0.973906528517171720077964012084452;
      x[2]=0.930157491355708226001207180059508;
      x[3]=0.865063366688984510732096688423493;
      x[4]=0.780817726586416897063717578345042;
      x[5]=0.679409568299024406234327365114874;
      x[6]=0.562757134668604683339000099272694;
      x[7]=0.433395394129247190799265943165784;
      x[8]=0.294392862701460198131126603103866;
      x[9]=0.148874338981631210884826001129720;
      x[10]=0.000000000000000000000000000000000;
      wkronrod[0]=0.011694638867371874278064396062192;
      wkronrod[1]=0.032558162307964727478818972459390;
      wkronrod[2]=0.054755896574351996031381300244580;
      wkronrod[3]=0.075039674810919952767043140916190;
      wkronrod[4]=0.093125454583697605535065465083366;
      wkronrod[5]=0.109387158802297641899210590325805;
      wkronrod[6]=0.123491976262065851077958109831074;
      wkronrod[7]=0.134709217311473325928054001771707;
      wkronrod[8]=0.142775938577060080797094273138717;
      wkronrod[9]=0.147739104901338491374841515972068;
      wkronrod[10]=0.149445554002916905664936468389821;
     }
//--- check
   if(n==31)
     {
      //--- change values
      ng=8;
      wgauss[0]=0.030753241996117268354628393577204;
      wgauss[1]=0.070366047488108124709267416450667;
      wgauss[2]=0.107159220467171935011869546685869;
      wgauss[3]=0.139570677926154314447804794511028;
      wgauss[4]=0.166269205816993933553200860481209;
      wgauss[5]=0.186161000015562211026800561866423;
      wgauss[6]=0.198431485327111576456118326443839;
      wgauss[7]=0.202578241925561272880620199967519;
      x[0]=0.998002298693397060285172840152271;
      x[1]=0.987992518020485428489565718586613;
      x[2]=0.967739075679139134257347978784337;
      x[3]=0.937273392400705904307758947710209;
      x[4]=0.897264532344081900882509656454496;
      x[5]=0.848206583410427216200648320774217;
      x[6]=0.790418501442465932967649294817947;
      x[7]=0.724417731360170047416186054613938;
      x[8]=0.650996741297416970533735895313275;
      x[9]=0.570972172608538847537226737253911;
      x[10]=0.485081863640239680693655740232351;
      x[11]=0.394151347077563369897207370981045;
      x[12]=0.299180007153168812166780024266389;
      x[13]=0.201194093997434522300628303394596;
      x[14]=0.101142066918717499027074231447392;
      x[15]=0.000000000000000000000000000000000;
      wkronrod[0]=0.005377479872923348987792051430128;
      wkronrod[1]=0.015007947329316122538374763075807;
      wkronrod[2]=0.025460847326715320186874001019653;
      wkronrod[3]=0.035346360791375846222037948478360;
      wkronrod[4]=0.044589751324764876608227299373280;
      wkronrod[5]=0.053481524690928087265343147239430;
      wkronrod[6]=0.062009567800670640285139230960803;
      wkronrod[7]=0.069854121318728258709520077099147;
      wkronrod[8]=0.076849680757720378894432777482659;
      wkronrod[9]=0.083080502823133021038289247286104;
      wkronrod[10]=0.088564443056211770647275443693774;
      wkronrod[11]=0.093126598170825321225486872747346;
      wkronrod[12]=0.096642726983623678505179907627589;
      wkronrod[13]=0.099173598721791959332393173484603;
      wkronrod[14]=0.100769845523875595044946662617570;
      wkronrod[15]=0.101330007014791549017374792767493;
     }
//--- check
   if(n==41)
     {
      //--- change values
      ng=10;
      wgauss[0]=0.017614007139152118311861962351853;
      wgauss[1]=0.040601429800386941331039952274932;
      wgauss[2]=0.062672048334109063569506535187042;
      wgauss[3]=0.083276741576704748724758143222046;
      wgauss[4]=0.101930119817240435036750135480350;
      wgauss[5]=0.118194531961518417312377377711382;
      wgauss[6]=0.131688638449176626898494499748163;
      wgauss[7]=0.142096109318382051329298325067165;
      wgauss[8]=0.149172986472603746787828737001969;
      wgauss[9]=0.152753387130725850698084331955098;
      x[0]=0.998859031588277663838315576545863;
      x[1]=0.993128599185094924786122388471320;
      x[2]=0.981507877450250259193342994720217;
      x[3]=0.963971927277913791267666131197277;
      x[4]=0.940822633831754753519982722212443;
      x[5]=0.912234428251325905867752441203298;
      x[6]=0.878276811252281976077442995113078;
      x[7]=0.839116971822218823394529061701521;
      x[8]=0.795041428837551198350638833272788;
      x[9]=0.746331906460150792614305070355642;
      x[10]=0.693237656334751384805490711845932;
      x[11]=0.636053680726515025452836696226286;
      x[12]=0.575140446819710315342946036586425;
      x[13]=0.510867001950827098004364050955251;
      x[14]=0.443593175238725103199992213492640;
      x[15]=0.373706088715419560672548177024927;
      x[16]=0.301627868114913004320555356858592;
      x[17]=0.227785851141645078080496195368575;
      x[18]=0.152605465240922675505220241022678;
      x[19]=0.076526521133497333754640409398838;
      x[20]=0.000000000000000000000000000000000;
      wkronrod[0]=0.003073583718520531501218293246031;
      wkronrod[1]=0.008600269855642942198661787950102;
      wkronrod[2]=0.014626169256971252983787960308868;
      wkronrod[3]=0.020388373461266523598010231432755;
      wkronrod[4]=0.025882133604951158834505067096153;
      wkronrod[5]=0.031287306777032798958543119323801;
      wkronrod[6]=0.036600169758200798030557240707211;
      wkronrod[7]=0.041668873327973686263788305936895;
      wkronrod[8]=0.046434821867497674720231880926108;
      wkronrod[9]=0.050944573923728691932707670050345;
      wkronrod[10]=0.055195105348285994744832372419777;
      wkronrod[11]=0.059111400880639572374967220648594;
      wkronrod[12]=0.062653237554781168025870122174255;
      wkronrod[13]=0.065834597133618422111563556969398;
      wkronrod[14]=0.068648672928521619345623411885368;
      wkronrod[15]=0.071054423553444068305790361723210;
      wkronrod[16]=0.073030690332786667495189417658913;
      wkronrod[17]=0.074582875400499188986581418362488;
      wkronrod[18]=0.075704497684556674659542775376617;
      wkronrod[19]=0.076377867672080736705502835038061;
      wkronrod[20]=0.076600711917999656445049901530102;
     }
//--- check
   if(n==51)
     {
      //--- change values
      ng=13;
      wgauss[0]=0.011393798501026287947902964113235;
      wgauss[1]=0.026354986615032137261901815295299;
      wgauss[2]=0.040939156701306312655623487711646;
      wgauss[3]=0.054904695975835191925936891540473;
      wgauss[4]=0.068038333812356917207187185656708;
      wgauss[5]=0.080140700335001018013234959669111;
      wgauss[6]=0.091028261982963649811497220702892;
      wgauss[7]=0.100535949067050644202206890392686;
      wgauss[8]=0.108519624474263653116093957050117;
      wgauss[9]=0.114858259145711648339325545869556;
      wgauss[10]=0.119455763535784772228178126512901;
      wgauss[11]=0.122242442990310041688959518945852;
      wgauss[12]=0.123176053726715451203902873079050;
      x[0]=0.999262104992609834193457486540341;
      x[1]=0.995556969790498097908784946893902;
      x[2]=0.988035794534077247637331014577406;
      x[3]=0.976663921459517511498315386479594;
      x[4]=0.961614986425842512418130033660167;
      x[5]=0.942974571228974339414011169658471;
      x[6]=0.920747115281701561746346084546331;
      x[7]=0.894991997878275368851042006782805;
      x[8]=0.865847065293275595448996969588340;
      x[9]=0.833442628760834001421021108693570;
      x[10]=0.797873797998500059410410904994307;
      x[11]=0.759259263037357630577282865204361;
      x[12]=0.717766406813084388186654079773298;
      x[13]=0.673566368473468364485120633247622;
      x[14]=0.626810099010317412788122681624518;
      x[15]=0.577662930241222967723689841612654;
      x[16]=0.526325284334719182599623778158010;
      x[17]=0.473002731445714960522182115009192;
      x[18]=0.417885382193037748851814394594572;
      x[19]=0.361172305809387837735821730127641;
      x[20]=0.303089538931107830167478909980339;
      x[21]=0.243866883720988432045190362797452;
      x[22]=0.183718939421048892015969888759528;
      x[23]=0.122864692610710396387359818808037;
      x[24]=0.061544483005685078886546392366797;
      x[25]=0.000000000000000000000000000000000;
      wkronrod[0]=0.001987383892330315926507851882843;
      wkronrod[1]=0.005561932135356713758040236901066;
      wkronrod[2]=0.009473973386174151607207710523655;
      wkronrod[3]=0.013236229195571674813656405846976;
      wkronrod[4]=0.016847817709128298231516667536336;
      wkronrod[5]=0.020435371145882835456568292235939;
      wkronrod[6]=0.024009945606953216220092489164881;
      wkronrod[7]=0.027475317587851737802948455517811;
      wkronrod[8]=0.030792300167387488891109020215229;
      wkronrod[9]=0.034002130274329337836748795229551;
      wkronrod[10]=0.037116271483415543560330625367620;
      wkronrod[11]=0.040083825504032382074839284467076;
      wkronrod[12]=0.042872845020170049476895792439495;
      wkronrod[13]=0.045502913049921788909870584752660;
      wkronrod[14]=0.047982537138836713906392255756915;
      wkronrod[15]=0.050277679080715671963325259433440;
      wkronrod[16]=0.052362885806407475864366712137873;
      wkronrod[17]=0.054251129888545490144543370459876;
      wkronrod[18]=0.055950811220412317308240686382747;
      wkronrod[19]=0.057437116361567832853582693939506;
      wkronrod[20]=0.058689680022394207961974175856788;
      wkronrod[21]=0.059720340324174059979099291932562;
      wkronrod[22]=0.060539455376045862945360267517565;
      wkronrod[23]=0.061128509717053048305859030416293;
      wkronrod[24]=0.061471189871425316661544131965264;
      wkronrod[25]=0.061580818067832935078759824240055;
     }
//--- check
   if(n==61)
     {
      //--- change values
      ng=15;
      wgauss[0]=0.007968192496166605615465883474674;
      wgauss[1]=0.018466468311090959142302131912047;
      wgauss[2]=0.028784707883323369349719179611292;
      wgauss[3]=0.038799192569627049596801936446348;
      wgauss[4]=0.048402672830594052902938140422808;
      wgauss[5]=0.057493156217619066481721689402056;
      wgauss[6]=0.065974229882180495128128515115962;
      wgauss[7]=0.073755974737705206268243850022191;
      wgauss[8]=0.080755895229420215354694938460530;
      wgauss[9]=0.086899787201082979802387530715126;
      wgauss[10]=0.092122522237786128717632707087619;
      wgauss[11]=0.096368737174644259639468626351810;
      wgauss[12]=0.099593420586795267062780282103569;
      wgauss[13]=0.101762389748405504596428952168554;
      wgauss[14]=0.102852652893558840341285636705415;
      x[0]=0.999484410050490637571325895705811;
      x[1]=0.996893484074649540271630050918695;
      x[2]=0.991630996870404594858628366109486;
      x[3]=0.983668123279747209970032581605663;
      x[4]=0.973116322501126268374693868423707;
      x[5]=0.960021864968307512216871025581798;
      x[6]=0.944374444748559979415831324037439;
      x[7]=0.926200047429274325879324277080474;
      x[8]=0.905573307699907798546522558925958;
      x[9]=0.882560535792052681543116462530226;
      x[10]=0.857205233546061098958658510658944;
      x[11]=0.829565762382768397442898119732502;
      x[12]=0.799727835821839083013668942322683;
      x[13]=0.767777432104826194917977340974503;
      x[14]=0.733790062453226804726171131369528;
      x[15]=0.697850494793315796932292388026640;
      x[16]=0.660061064126626961370053668149271;
      x[17]=0.620526182989242861140477556431189;
      x[18]=0.579345235826361691756024932172540;
      x[19]=0.536624148142019899264169793311073;
      x[20]=0.492480467861778574993693061207709;
      x[21]=0.447033769538089176780609900322854;
      x[22]=0.400401254830394392535476211542661;
      x[23]=0.352704725530878113471037207089374;
      x[24]=0.304073202273625077372677107199257;
      x[25]=0.254636926167889846439805129817805;
      x[26]=0.204525116682309891438957671002025;
      x[27]=0.153869913608583546963794672743256;
      x[28]=0.102806937966737030147096751318001;
      x[29]=0.051471842555317695833025213166723;
      x[30]=0.000000000000000000000000000000000;
      wkronrod[0]=0.001389013698677007624551591226760;
      wkronrod[1]=0.003890461127099884051267201844516;
      wkronrod[2]=0.006630703915931292173319826369750;
      wkronrod[3]=0.009273279659517763428441146892024;
      wkronrod[4]=0.011823015253496341742232898853251;
      wkronrod[5]=0.014369729507045804812451432443580;
      wkronrod[6]=0.016920889189053272627572289420322;
      wkronrod[7]=0.019414141193942381173408951050128;
      wkronrod[8]=0.021828035821609192297167485738339;
      wkronrod[9]=0.024191162078080601365686370725232;
      wkronrod[10]=0.026509954882333101610601709335075;
      wkronrod[11]=0.028754048765041292843978785354334;
      wkronrod[12]=0.030907257562387762472884252943092;
      wkronrod[13]=0.032981447057483726031814191016854;
      wkronrod[14]=0.034979338028060024137499670731468;
      wkronrod[15]=0.036882364651821229223911065617136;
      wkronrod[16]=0.038678945624727592950348651532281;
      wkronrod[17]=0.040374538951535959111995279752468;
      wkronrod[18]=0.041969810215164246147147541285970;
      wkronrod[19]=0.043452539701356069316831728117073;
      wkronrod[20]=0.044814800133162663192355551616723;
      wkronrod[21]=0.046059238271006988116271735559374;
      wkronrod[22]=0.047185546569299153945261478181099;
      wkronrod[23]=0.048185861757087129140779492298305;
      wkronrod[24]=0.049055434555029778887528165367238;
      wkronrod[25]=0.049795683427074206357811569379942;
      wkronrod[26]=0.050405921402782346840893085653585;
      wkronrod[27]=0.050881795898749606492297473049805;
      wkronrod[28]=0.051221547849258772170656282604944;
      wkronrod[29]=0.051426128537459025933862879215781;
      wkronrod[30]=0.051494729429451567558340433647099;
     }
//--- copy nodes
   for(i=n-1; i>=n/2; i--)
      x[i]=-x[n-1-i];
//--- copy Kronrod weights
   for(i=n-1; i>=n/2; i--)
      wkronrod[i]=wkronrod[n-1-i];
//--- copy Gauss weights
   for(i=ng-1; i>=0; i--)
     {
      wgauss[n-2-2*i]=wgauss[i];
      wgauss[1+2*i]=wgauss[i];
     }
   for(i=0; i<=n/2; i++)
      wgauss[2*i]=0;
//--- reorder
   CTSort::TagSort(x,n,p1,p2);
   for(i=0; i<=n-1; i++)
     {
      //--- swap
      tmp=wkronrod[i];
      wkronrod[i]=wkronrod[p2[i]];
      wkronrod[p2[i]]=tmp;
      tmp=wgauss[i];
      wgauss[i]=wgauss[p2[i]];
      wgauss[p2[i]]=tmp;
     }
  }
//+------------------------------------------------------------------+
//| Integration report:                                              |
//| * TerminationType = completetion code:                           |
//|     * -5    non-convergence of Gauss-Kronrod nodes               |
//|             calculation subroutine.                              |
//|     * -1    incorrect parameters were specified                  |
//|     *  1    OK                                                   |
//| * Rep.NFEV countains number of function calculations             |
//| * Rep.NIntervals contains number of intervals [a,b]              |
//|   was partitioned into.                                          |
//+------------------------------------------------------------------+
class CAutoGKReport
  {
public:
   //--- variables
   int               m_terminationtype;
   int               m_nfev;
   int               m_nintervals;
   //--- constructor, destructor
                     CAutoGKReport(void) { ZeroMemory(this); }
                    ~CAutoGKReport(void) {}
   //--- copy
   void              Copy(CAutoGKReport &obj);
  };
//+------------------------------------------------------------------+
//| Copy                                                             |
//+------------------------------------------------------------------+
void CAutoGKReport::Copy(CAutoGKReport &obj)
  {
//--- copy variables
   m_terminationtype=obj.m_terminationtype;
   m_nfev=obj.m_nfev;
   m_nintervals=obj.m_nintervals;
  }
//+------------------------------------------------------------------+
//| Integration report:                                              |
//| * TerminationType = completetion code:                           |
//|     * -5    non-convergence of Gauss-Kronrod nodes               |
//|             calculation subroutine.                              |
//|     * -1    incorrect parameters were specified                  |
//|     *  1    OK                                                   |
//| * Rep.NFEV countains number of function calculations             |
//| * Rep.NIntervals contains number of intervals [a,b]              |
//|   was partitioned into.                                          |
//+------------------------------------------------------------------+
class CAutoGKReportShell
  {
private:
   CAutoGKReport     m_innerobj;

public:
   //--- constructors, destructor
                     CAutoGKReportShell(void) {}
                     CAutoGKReportShell(CAutoGKReport &obj) { m_innerobj.Copy(obj); }
                    ~CAutoGKReportShell(void) {}
   //--- methods
   int               GetTerminationType(void);
   void              SetTerminationType(const int i);
   int               GetNFev(void);
   void              SetNFev(const int i);
   int               GetNIntervals(void);
   void              SetNIntervals(const int i);
   CAutoGKReport    *GetInnerObj(void);
  };
//+------------------------------------------------------------------+
//| Returns the value of the variable nfev                           |
//+------------------------------------------------------------------+
int CAutoGKReportShell::GetNFev(void)
  {
   return(m_innerobj.m_nfev);
  }
//+------------------------------------------------------------------+
//| Changing the value of the variable nfev                          |
//+------------------------------------------------------------------+
void CAutoGKReportShell::SetNFev(const int i)
  {
   m_innerobj.m_nfev=i;
  }
//+------------------------------------------------------------------+
//| Returns the value of the variable terminationtype                |
//+------------------------------------------------------------------+
int CAutoGKReportShell::GetTerminationType(void)
  {
   return(m_innerobj.m_terminationtype);
  }
//+------------------------------------------------------------------+
//| Changing the value of the variable terminationtype               |
//+------------------------------------------------------------------+
void CAutoGKReportShell::SetTerminationType(const int i)
  {
   m_innerobj.m_terminationtype=i;
  }
//+------------------------------------------------------------------+
//| Returns the value of the variable nintervals                     |
//+------------------------------------------------------------------+
int CAutoGKReportShell::GetNIntervals(void)
  {
   return(m_innerobj.m_nintervals);
  }
//+------------------------------------------------------------------+
//| Changing the value of the variable nintervals                    |
//+------------------------------------------------------------------+
void CAutoGKReportShell::SetNIntervals(const int i)
  {
   m_innerobj.m_nintervals=i;
  }
//+------------------------------------------------------------------+
//| Return object of class                                           |
//+------------------------------------------------------------------+
CAutoGKReport *CAutoGKReportShell::GetInnerObj(void)
  {
   return(GetPointer(m_innerobj));
  }
//+------------------------------------------------------------------+
//| Auxiliary class for CAutoGK                                      |
//+------------------------------------------------------------------+
struct CAutoGKInternalState
  {
   //--- variables
   double            m_a;
   double            m_b;
   double            m_eps;
   double            m_xwidth;
   double            m_x;
   double            m_f;
   int               m_info;
   double            m_r;
   int               m_heapsize;
   int               m_heapwidth;
   int               m_heapused;
   double            m_sumerr;
   double            m_sumabs;
   int               m_n;
   RCommState        m_rstate;
   //--- arrays
   double            m_qn[];
   double            m_wg[];
   double            m_wk[];
   double            m_wr[];
   //--- matrix
   CMatrixDouble     m_heap;
   //--- constructor, destructor
                     CAutoGKInternalState(void);
                    ~CAutoGKInternalState(void) {}
   //--- copy
   void              Copy(CAutoGKInternalState &obj);
  };
//+------------------------------------------------------------------+
//| Constructor                                                      |
//+------------------------------------------------------------------+
CAutoGKInternalState::CAutoGKInternalState(void)
  {
   m_a=0;
   m_b=0;
   m_eps=0;
   m_xwidth=0;
   m_x=0;
   m_f=0;
   m_info=0;
   m_r=0;
   m_heapsize=0;
   m_heapwidth=0;
   m_heapused=0;
   m_sumerr=0;
   m_sumabs=0;
   m_n=0;
  }
//+------------------------------------------------------------------+
//| Copy                                                             |
//+------------------------------------------------------------------+
void CAutoGKInternalState::Copy(CAutoGKInternalState &obj)
  {
//--- copy variables
   m_a=obj.m_a;
   m_b=obj.m_b;
   m_eps=obj.m_eps;
   m_xwidth=obj.m_xwidth;
   m_x=obj.m_x;
   m_f=obj.m_f;
   m_info=obj.m_info;
   m_r=obj.m_r;
   m_heapsize=obj.m_heapsize;
   m_heapwidth=obj.m_heapwidth;
   m_heapused=obj.m_heapused;
   m_sumerr=obj.m_sumerr;
   m_sumabs=obj.m_sumabs;
   m_n=obj.m_n;
   m_rstate.Copy(obj.m_rstate);
//--- copy arrays
   ArrayCopy(m_qn,obj.m_qn);
   ArrayCopy(m_wg,obj.m_wg);
   ArrayCopy(m_wk,obj.m_wk);
   ArrayCopy(m_wr,obj.m_wr);
//--- copy matrix
   m_heap=obj.m_heap;
  }
//+------------------------------------------------------------------+
//| This structure stores state of the integration algorithm.        |
//| Although this class has public fields,  they are not intended for|
//| external use. You should use ALGLIB functions to work with this  |
//| class:                                                           |
//| * autogksmooth()/AutoGKSmoothW()/... to create objects           |
//| * autogkintegrate() to begin integration                         |
//| * autogkresults() to get results                                 |
//+------------------------------------------------------------------+
class CAutoGKState
  {
public:
   //--- variables
   double            m_a;
   double            m_b;
   double            m_alpha;
   double            m_beta;
   double            m_xwidth;
   double            m_x;
   double            m_xminusa;
   double            m_bminusx;
   bool              m_needf;
   double            m_f;
   int               m_wrappermode;
   double            m_v;
   int               m_terminationtype;
   int               m_nfev;
   int               m_nintervals;
   CAutoGKInternalState m_internalstate;
   RCommState        m_rstate;
   //--- constructor, destructor
                     CAutoGKState(void);
                    ~CAutoGKState(void) {}
   //--- copy
   void              Copy(CAutoGKState &obj);
  };
//+------------------------------------------------------------------+
//| Constructor                                                      |
//+------------------------------------------------------------------+
CAutoGKState::CAutoGKState(void)
  {
   m_a=0;
   m_b=0;
   m_alpha=0;
   m_beta=0;
   m_xwidth=0;
   m_x=0;
   m_xminusa=0;
   m_bminusx=0;
   m_needf=false;
   m_f=0;
   m_wrappermode=0;
   m_v=0;
   m_terminationtype=0;
   m_nfev=0;
   m_nintervals=0;
  }
//+------------------------------------------------------------------+
//| Copy                                                             |
//+------------------------------------------------------------------+
void CAutoGKState::Copy(CAutoGKState &obj)
  {
//--- copy variables
   m_a=obj.m_a;
   m_b=obj.m_b;
   m_alpha=obj.m_alpha;
   m_beta=obj.m_beta;
   m_xwidth=obj.m_xwidth;
   m_x=obj.m_x;
   m_xminusa=obj.m_xminusa;
   m_bminusx=obj.m_bminusx;
   m_needf=obj.m_needf;
   m_f=obj.m_f;
   m_wrappermode=obj.m_wrappermode;
   m_v=obj.m_v;
   m_terminationtype=obj.m_terminationtype;
   m_nfev=obj.m_nfev;
   m_nintervals=obj.m_nintervals;
   m_internalstate.Copy(obj.m_internalstate);
   m_rstate.Copy(obj.m_rstate);
  }
//+------------------------------------------------------------------+
//| This structure stores state of the integration algorithm.        |
//| Although this class has public fields,  they are not intended for|
//| external use. You should use ALGLIB functions to work with this  |
//| class:                                                           |
//| * autogksmooth()/AutoGKSmoothW()/... to create objects           |
//| * autogkintegrate() to begin integration                         |
//| * autogkresults() to get results                                 |
//+------------------------------------------------------------------+
class CAutoGKStateShell
  {
private:
   CAutoGKState      m_innerobj;

public:
   //--- constructors, destructor
                     CAutoGKStateShell(void) {}
                     CAutoGKStateShell(CAutoGKState &obj) { m_innerobj.Copy(obj); }
                    ~CAutoGKStateShell(void) {}
   //--- methods
   bool              GetNeedF(void);
   void              SetNeedF(const bool b);
   double            GetX(void);
   void              SetX(const double d);
   double            GetXMinusA(void);
   void              SetXMinusA(const double d);
   double            GetBMinusX(void);
   void              SetBMinusX(const double d);
   double            GetF(void);
   void              SetF(const double d);
   CAutoGKState     *GetInnerObj(void);
  };
//+------------------------------------------------------------------+
//| Returns the value of the variable needf                          |
//+------------------------------------------------------------------+
bool CAutoGKStateShell::GetNeedF(void)
  {
   return(m_innerobj.m_needf);
  }
//+------------------------------------------------------------------+
//| Changing the value of the variable needf                         |
//+------------------------------------------------------------------+
void CAutoGKStateShell::SetNeedF(const bool b)
  {
   m_innerobj.m_needf=b;
  }
//+------------------------------------------------------------------+
//| Returns the value of the variable x                              |
//+------------------------------------------------------------------+
double CAutoGKStateShell::GetX(void)
  {
   return(m_innerobj.m_x);
  }
//+------------------------------------------------------------------+
//| Changing the value of the variable x                             |
//+------------------------------------------------------------------+
void CAutoGKStateShell::SetX(const double d)
  {
   m_innerobj.m_x=d;
  }
//+------------------------------------------------------------------+
//| Returns the value of the variable xminusa                        |
//+------------------------------------------------------------------+
double CAutoGKStateShell::GetXMinusA(void)
  {
   return(m_innerobj.m_xminusa);
  }
//+------------------------------------------------------------------+
//| Changing the value of the variable xminusa                       |
//+------------------------------------------------------------------+
void CAutoGKStateShell::SetXMinusA(const double d)
  {
   m_innerobj.m_xminusa=d;
  }
//+------------------------------------------------------------------+
//| Returns the value of the variable bminusx                        |
//+------------------------------------------------------------------+
double CAutoGKStateShell::GetBMinusX(void)
  {
   return(m_innerobj.m_bminusx);
  }
//+------------------------------------------------------------------+
//| Changing the value of the variable bminusx                       |
//+------------------------------------------------------------------+
void CAutoGKStateShell::SetBMinusX(const double d)
  {
   m_innerobj.m_bminusx=d;
  }
//+------------------------------------------------------------------+
//| Returns the value of the variable f                              |
//+------------------------------------------------------------------+
double CAutoGKStateShell::GetF(void)
  {
   return(m_innerobj.m_f);
  }
//+------------------------------------------------------------------+
//| Changing the value of the variable f                             |
//+------------------------------------------------------------------+
void CAutoGKStateShell::SetF(const double d)
  {
   m_innerobj.m_f=d;
  }
//+------------------------------------------------------------------+
//| Return object of class                                           |
//+------------------------------------------------------------------+
CAutoGKState *CAutoGKStateShell::GetInnerObj(void)
  {
   return(GetPointer(m_innerobj));
  }
//+------------------------------------------------------------------+
//| Auto Gauss-Kronrod                                               |
//+------------------------------------------------------------------+
class CAutoGK
  {
public:
   //--- class constant
   static const int  m_maxsubintervals;
   //--- public methods
   static void       AutoGKSmooth(const double a,const double b,CAutoGKState &state);
   static void       AutoGKSmoothW(const double a,const double b,const double xwidth,CAutoGKState &state);
   static void       AutoGKSingular(const double a,const double b,const double alpha,const double beta,CAutoGKState &state);
   static void       AutoGKResults(CAutoGKState &state,double &v,CAutoGKReport &rep);
   static bool       AutoGKIteration(CAutoGKState &state);

private:
   //--- private methods
   static void       AutoGKInternalPrepare(const double a,const double b,const double eps,const double xwidth,CAutoGKInternalState &state);
   static void       MHeapPop(CMatrixDouble &heap,const int heapsize,const int heapwidth);
   static void       MHeapPush(CMatrixDouble &heap,const int heapsize,const int heapwidth);
   static void       MHeapResize(CMatrixDouble &heap,int &heapsize,const int newheapsize,const int heapwidth);
   static bool       AutoGKInternalIteration(CAutoGKInternalState &state);
   //--- auxiliary functions for AutoGKInternalIteration
   static void       Func_Internal_lbl_rcomm(CAutoGKInternalState &state,int i,int j,int ns,int info,double c1,double c2,double intg,double intk,double inta,double v,double ta,double tb,double qeps);
   static bool       Func_Internal_lbl_5(CAutoGKInternalState &state,int &i,int &j,int &ns,int &info,double &c1,double &c2,double &intg,double &intk,double &inta,double &v,double &ta,double &tb,double &qeps);
   static bool       Func_Internal_lbl_7(CAutoGKInternalState &state,int &i,int &j,int &ns,int &info,double &c1,double &c2,double &intg,double &intk,double &inta,double &v,double &ta,double &tb,double &qeps);
   static bool       Func_Internal_lbl_8(CAutoGKInternalState &state,int &i,int &j,int &ns,int &info,double &c1,double &c2,double &intg,double &intk,double &inta,double &v,double &ta,double &tb,double &qeps);
   static bool       Func_Internal_lbl_11(CAutoGKInternalState &state,int &i,int &j,int &ns,int &info,double &c1,double &c2,double &intg,double &intk,double &inta,double &v,double &ta,double &tb,double &qeps);
   static bool       Func_Internal_lbl_14(CAutoGKInternalState &state,int &i,int &j,int &ns,int &info,double &c1,double &c2,double &intg,double &intk,double &inta,double &v,double &ta,double &tb,double &qeps);
   static bool       Func_Internal_lbl_16(CAutoGKInternalState &state,int &i,int &j,int &ns,int &info,double &c1,double &c2,double &intg,double &intk,double &inta,double &v,double &ta,double &tb,double &qeps);
   static bool       Func_Internal_lbl_19(CAutoGKInternalState &state,int &i,int &j,int &ns,int &info,double &c1,double &c2,double &intg,double &intk,double &inta,double &v,double &ta,double &tb,double &qeps);
   //--- auxiliary functions for AutoGKIteration
   static void       Func_lbl_rcomm(CAutoGKState &state,double s,double tmp,double eps,double a,double b,double x,double t,double alpha,double beta,double v1,double v2);
   static bool       Func_lbl_5(CAutoGKState &state,double &s,double &tmp,double &eps,double &a,double &b,double &x,double &t,double &alpha,double &beta,double &v1,double &v2);
   static bool       Func_lbl_9(CAutoGKState &state,double &s,double &tmp,double &eps,double &a,double &b,double &x,double &t,double &alpha,double &beta,double &v1,double &v2);
   static bool       Func_lbl_11(CAutoGKState &state,double &s,double &tmp,double &eps,double &a,double &b,double &x,double &t,double &alpha,double &beta,double &v1,double &v2);
  };
//+------------------------------------------------------------------+
//| Initialize constant                                              |
//+------------------------------------------------------------------+
const int CAutoGK::m_maxsubintervals=10000;
//+------------------------------------------------------------------+
//| Integration of a smooth function F(x) on a finite interval [a,b].|
//| Fast-convergent algorithm based on a Gauss-Kronrod formula is    |
//| used. Result is calculated with accuracy close to the machine    |
//| precision.                                                       |
//| Algorithm works well only with smooth integrands. It may be used |
//| with continuous non-smooth integrands, but with less performance.|
//| It should never be used with integrands which have integrable    |
//| singularities at lower or upper limits - algorithm may crash.    |
//| Use AutoGKSingular in such cases.                                |
//| INPUT PARAMETERS:                                                |
//|     A, B    -   interval boundaries (A<B, A=B or A>B)            |
//| OUTPUT PARAMETERS                                                |
//|     State   -   structure which stores algorithm state           |
//| SEE ALSO                                                         |
//|     AutoGKSmoothW, AutoGKSingular, AutoGKResults.                |
//+------------------------------------------------------------------+
void CAutoGK::AutoGKSmooth(const double a,const double b,
                           CAutoGKState &state)
  {
//--- check
   if(!CAp::Assert(CMath::IsFinite(a),__FUNCTION__+": A is not finite!"))
      return;
//--- check
   if(!CAp::Assert(CMath::IsFinite(b),__FUNCTION__+": B is not finite!"))
      return;
//--- function call
   AutoGKSmoothW(a,b,0.0,state);
  }
//+------------------------------------------------------------------+
//| Integration of a smooth function F(x) on a finite interval [a,b].|
//| This subroutine is same as AutoGKSmooth(), but it guarantees that|
//| interval [a,b] is partitioned into subintervals which have width |
//| at most XWidth.                                                  |
//| Subroutine can be used when integrating nearly-constant function |
//| with narrow "bumps" (about XWidth wide). If "bumps" are too      |
//| narrow, AutoGKSmooth subroutine can overlook them.               |
//| INPUT PARAMETERS:                                                |
//|     A, B    -   interval boundaries (A<B, A=B or A>B)            |
//| OUTPUT PARAMETERS                                                |
//|     State   -   structure which stores algorithm state           |
//| SEE ALSO                                                         |
//|     AutoGKSmooth, AutoGKSingular, AutoGKResults.                 |
//+------------------------------------------------------------------+
void CAutoGK::AutoGKSmoothW(const double a,const double b,
                            const double xwidth,CAutoGKState &state)
  {
//--- check
   if(!CAp::Assert(CMath::IsFinite(a),__FUNCTION__+": A is not finite!"))
      return;
//--- check
   if(!CAp::Assert(CMath::IsFinite(b),__FUNCTION__+": B is not finite!"))
      return;
//--- check
   if(!CAp::Assert(CMath::IsFinite(xwidth),__FUNCTION__+": XWidth is not finite!"))
      return;
//--- initialization
   state.m_wrappermode=0;
   state.m_a=a;
   state.m_b=b;
   state.m_xwidth=xwidth;
   state.m_needf=false;
//--- allocation
   state.m_rstate.ra.Resize(11);
   state.m_rstate.stage=-1;
  }
//+------------------------------------------------------------------+
//| Integration on a finite interval [A,B].                          |
//| Integrand have integrable singularities at A/B.                  |
//| F(X) must diverge as "(x-A)^alpha" at A, as "(B-x)^beta" at B,   |
//| with known alpha/beta (alpha>-1, beta>-1). If alpha/beta are not |
//| known, estimates from below can be used (but these estimates     |
//| should be greater than -1 too).                                  |
//| One of alpha/beta variables (or even both alpha/beta) may be     |
//| equal to 0, which means than function F(x) is non-singular at    |
//| A/B. Anyway (singular at bounds or not), function F(x) is        |
//| supposed to be continuous on (A,B).                              |
//| Fast-convergent algorithm based on a Gauss-Kronrod formula is    |
//| used. Result is calculated with accuracy close to the machine    |
//| precision.                                                       |
//| INPUT PARAMETERS:                                                |
//|     A, B    -   interval boundaries (A<B, A=B or A>B)            |
//|     Alpha   -   power-law coefficient of the F(x) at A,          |
//|                 Alpha>-1                                         |
//|     Beta    -   power-law coefficient of the F(x) at B,          |
//|                 Beta>-1                                          |
//| OUTPUT PARAMETERS                                                |
//|     State   -   structure which stores algorithm state           |
//| SEE ALSO                                                         |
//|     AutoGKSmooth, AutoGKSmoothW, AutoGKResults.                  |
//+------------------------------------------------------------------+
void CAutoGK::AutoGKSingular(const double a,const double b,
                             const double alpha,const double beta,
                             CAutoGKState &state)
  {
//--- check
   if(!CAp::Assert(CMath::IsFinite(a),__FUNCTION__+": A is not finite!"))
      return;
//--- check
   if(!CAp::Assert(CMath::IsFinite(b),__FUNCTION__+": B is not finite!"))
      return;
//--- check
   if(!CAp::Assert(CMath::IsFinite(alpha),__FUNCTION__+": Alpha is not finite!"))
      return;
//--- check
   if(!CAp::Assert(CMath::IsFinite(beta),__FUNCTION__+": Beta is not finite!"))
      return;
//--- initialization
   state.m_wrappermode=1;
   state.m_a=a;
   state.m_b=b;
   state.m_alpha=alpha;
   state.m_beta=beta;
   state.m_xwidth=0.0;
   state.m_needf=false;
//--- allocation
   state.m_rstate.ra.Resize(11);
   state.m_rstate.stage=-1;
  }
//+------------------------------------------------------------------+
//| Adaptive integration results                                     |
//| Called after AutoGKIteration returned False.                     |
//| Input parameters:                                                |
//|     State   -   algorithm state (used by AutoGKIteration).       |
//| Output parameters:                                               |
//|     V       -   integral(f(x)dx,a,b)                             |
//|     Rep     -   optimization report (see AutoGKReport            |
//|                 description)                                     |
//+------------------------------------------------------------------+
void CAutoGK::AutoGKResults(CAutoGKState &state,double &v,CAutoGKReport &rep)
  {
//--- change values
   v=state.m_v;
   rep.m_terminationtype=state.m_terminationtype;
   rep.m_nfev=state.m_nfev;
   rep.m_nintervals=state.m_nintervals;
  }
//+------------------------------------------------------------------+
//| Internal AutoGK subroutine                                       |
//| eps<0   - error                                                  |
//| eps=0   - automatic eps selection                                |
//| width<0 -   error                                                |
//| width=0 -   no width requirements                                |
//+------------------------------------------------------------------+
void CAutoGK::AutoGKInternalPrepare(const double a,const double b,
                                    const double eps,const double xwidth,
                                    CAutoGKInternalState &state)
  {
//--- Save settings
   state.m_a=a;
   state.m_b=b;
   state.m_eps=eps;
   state.m_xwidth=xwidth;
//--- Prepare RComm structure
   state.m_rstate.ia.Resize(4);
   state.m_rstate.ra.Resize(9);
   state.m_rstate.stage=-1;
  }
//+------------------------------------------------------------------+
//| Internal subroutine                                              |
//+------------------------------------------------------------------+
void CAutoGK::MHeapPop(CMatrixDouble &heap,const int heapsize,
                       const int heapwidth)
  {
//--- create variables
   int    i=0;
   int    p=0;
   double t=0;
   int    maxcp=0;
//--- check
   if(heapsize==1)
      return;
//--- initialization
   for(i=0; i<=heapwidth-1; i++)
     {
      t=heap[heapsize-1][i];
      heap.Set(heapsize-1,i,heap[0][i]);
      heap.Set(0,i,t);
     }
   p=0;
//--- cycle
   while(2*p+1<heapsize-1)
     {
      maxcp=2*p+1;
      //--- check
      if(2*p+2<heapsize-1)
        {
         //--- check
         if(heap[2*p+2][0]>heap[2*p+1][0])
            maxcp=2*p+2;
        }
      //--- check
      if(heap[p][0]<heap[maxcp][0])
        {
         for(i=0; i<=heapwidth-1; i++)
           {
            //--- change values
            t=heap[p][i];
            heap.Set(p,i,heap[maxcp][i]);
            heap.Set(maxcp,i,t);
           }
         p=maxcp;
        }
      else
         break;
     }
  }
//+------------------------------------------------------------------+
//| Internal subroutine                                              |
//+------------------------------------------------------------------+
void CAutoGK::MHeapPush(CMatrixDouble &heap,const int heapsize,
                        const int heapwidth)
  {
//--- create variables
   int    i=0;
   int    p=0;
   double t=0;
   int    parent=0;
//--- check
   if(heapsize==0)
      return;
   p=heapsize;
//--- cycle
   while(p!=0)
     {
      parent=(p-1)/2;
      //--- check
      if(heap[p][0]>heap[parent][0])
        {
         for(i=0; i<=heapwidth-1; i++)
           {
            //--- change values
            t=heap[p][i];
            heap.Set(p,i,heap[parent][i]);
            heap.Set(parent,i,t);
           }
         p=parent;
        }
      else
         break;
     }
  }
//+------------------------------------------------------------------+
//| Internal subroutine                                              |
//+------------------------------------------------------------------+
void CAutoGK::MHeapResize(CMatrixDouble &heap,int &heapsize,
                          const int newheapsize,const int heapwidth)
  {
//--- create variables
   int i=0;
   int i_=0;
//--- create matrix
   CMatrixDouble tmp;
//--- allocation
   tmp.Resize(heapsize,heapwidth);
//--- copy
   for(i=0; i<=heapsize-1; i++)
     {
      for(i_=0; i_<=heapwidth-1; i_++)
         tmp.Set(i,i_,heap[i][i_]);
     }
//--- allocation
   heap.Resize(newheapsize,heapwidth);
//--- copy
   for(i=0; i<=heapsize-1; i++)
     {
      for(i_=0; i_<=heapwidth-1; i_++)
         heap.Set(i,i_,tmp[i][i_]);
     }
//--- change value
   heapsize=newheapsize;
  }
//+------------------------------------------------------------------+
//| Internal AutoGK subroutine                                       |
//+------------------------------------------------------------------+
bool CAutoGK::AutoGKInternalIteration(CAutoGKInternalState &state)
  {
//--- create variables
   double c1=0;
   double c2=0;
   int    i=0;
   int    j=0;
   double intg=0;
   double intk=0;
   double inta=0;
   double v=0;
   double ta=0;
   double tb=0;
   int    ns=0;
   double qeps=0;
   int    info=0;
//--- This code initializes locals by:
//--- * random values determined during code
//---   generation - on first subroutine call
//--- * values from previous call - on subsequent calls
   if(state.m_rstate.stage>=0)
     {
      //--- initialization
      i=state.m_rstate.ia[0];
      j=state.m_rstate.ia[1];
      ns=state.m_rstate.ia[2];
      info=state.m_rstate.ia[3];
      c1=state.m_rstate.ra[0];
      c2=state.m_rstate.ra[1];
      intg=state.m_rstate.ra[2];
      intk=state.m_rstate.ra[3];
      inta=state.m_rstate.ra[4];
      v=state.m_rstate.ra[5];
      ta=state.m_rstate.ra[6];
      tb=state.m_rstate.ra[7];
      qeps=state.m_rstate.ra[8];
     }
   else
     {
      //--- initialization
      i=497;
      j=-271;
      ns=-581;
      info=745;
      c1=-533;
      c2=-77;
      intg=678;
      intk=-293;
      inta=316;
      v=647;
      ta=-756;
      tb=830;
      qeps=-871;
     }
//--- check
   if(state.m_rstate.stage==0)
     {
      //--- change value
      v=state.m_f;
      //--- Gauss-Kronrod formula
      intk=intk+v*state.m_wk[i];
      //--- check
      if(i%2==1)
         intg=intg+v*state.m_wg[i];
      //--- Integral |F(x)|
      //--- Use rectangles method
      inta=inta+MathAbs(v)*state.m_wr[i];
      i=i+1;
      //--- function call, return result
      return(Func_Internal_lbl_5(state,i,j,ns,info,c1,c2,intg,intk,inta,v,ta,tb,qeps));
     }
//--- check
   if(state.m_rstate.stage==1)
     {
      //--- change value
      v=state.m_f;
      //--- Gauss-Kronrod formula
      intk=intk+v*state.m_wk[i];
      //--- check
      if(i%2==1)
         intg=intg+v*state.m_wg[i];
      //--- Integral |F(x)|
      //--- Use rectangles method
      inta=inta+MathAbs(v)*state.m_wr[i];
      i=i+1;
      //--- function call, return result
      return(Func_Internal_lbl_11(state,i,j,ns,info,c1,c2,intg,intk,inta,v,ta,tb,qeps));
     }
//--- check
   if(state.m_rstate.stage==2)
     {
      //--- change value
      v=state.m_f;
      //--- Gauss-Kronrod formula
      intk=intk+v*state.m_wk[i];
      //--- check
      if(i%2==1)
         intg=intg+v*state.m_wg[i];
      //--- Integral |F(x)|
      //--- Use rectangles method
      inta=inta+MathAbs(v)*state.m_wr[i];
      i=i+1;
      //--- function call, return result
      return(Func_Internal_lbl_19(state,i,j,ns,info,c1,c2,intg,intk,inta,v,ta,tb,qeps));
     }
//--- Routine body
//--- initialize quadratures.
//--- use 15-point Gauss-Kronrod formula.
   state.m_n=15;
   CGaussKronrodQ::GKQGenerateGaussLegendre(state.m_n,info,state.m_qn,state.m_wk,state.m_wg);
//--- check
   if(info<0)
     {
      //--- change values
      state.m_info=-5;
      state.m_r=0;
      //--- return result
      return(false);
     }
//--- allocation
   ArrayResize(state.m_wr,state.m_n);
   for(i=0; i<=state.m_n-1; i++)
     {
      //--- check
      if(i==0)
        {
         state.m_wr[i]=0.5*MathAbs(state.m_qn[1]-state.m_qn[0]);
         continue;
        }
      //--- check
      if(i==state.m_n-1)
        {
         state.m_wr[state.m_n-1]=0.5*MathAbs(state.m_qn[state.m_n-1]-state.m_qn[state.m_n-2]);
         continue;
        }
      state.m_wr[i]=0.5*MathAbs(state.m_qn[i-1]-state.m_qn[i+1]);
     }
//--- special case
   if(state.m_a==state.m_b)
     {
      //--- change values
      state.m_info=1;
      state.m_r=0;
      //--- return result
      return(false);
     }
//--- test parameters
   if(state.m_eps<0.0 || state.m_xwidth<0.0)
     {
      //--- change values
      state.m_info=-1;
      state.m_r=0;
      //--- return result
      return(false);
     }
   state.m_info=1;
//--- check
   if(state.m_eps==0.0)
      state.m_eps=100000*CMath::m_machineepsilon;
//--- First,prepare heap
//--- * column 0   -   absolute error
//--- * column 1   -   integral of a F(x) (calculated using Kronrod extension nodes)
//--- * column 2   -   integral of a |F(x)| (calculated using modified rect. method)
//--- * column 3   -   left boundary of a subinterval
//--- * column 4   -   right boundary of a subinterval
   if(state.m_xwidth!=0.0)
     {
      //--- maximum subinterval should be no more than XWidth.
      //--- so we create Ceil((B-A)/XWidth)+1 small subintervals
      ns=(int)MathCeil(MathAbs(state.m_b-state.m_a)/state.m_xwidth)+1;
      state.m_heapsize=ns;
      state.m_heapused=ns;
      state.m_heapwidth=5;
      state.m_heap.Resize(state.m_heapsize,state.m_heapwidth);
      state.m_sumerr=0;
      state.m_sumabs=0;
      j=0;
      //--- function call, return result
      return(Func_Internal_lbl_8(state,i,j,ns,info,c1,c2,intg,intk,inta,v,ta,tb,qeps));
     }
//--- no maximum width requirements
//--- start from one big subinterval
   state.m_heapwidth=5;
   state.m_heapsize=1;
   state.m_heapused=1;
   state.m_heap.Resize(state.m_heapsize,state.m_heapwidth);
   c1=0.5*(state.m_b-state.m_a);
   c2=0.5*(state.m_b+state.m_a);
   intg=0;
   intk=0;
   inta=0;
   i=0;
//--- function call, return result
   return(Func_Internal_lbl_5(state,i,j,ns,info,c1,c2,intg,intk,inta,v,ta,tb,qeps));
  }
//+------------------------------------------------------------------+
//| Auxiliary function for AutoGKInternalIteration. Is a product to  |
//| get rid of the operator unconditional jump goto.                 |
//+------------------------------------------------------------------+
void CAutoGK::Func_Internal_lbl_rcomm(CAutoGKInternalState &state,
                                      int i,int j,int ns,int info,
                                      double c1,double c2,double intg,
                                      double intk,double inta,double v,
                                      double ta,double tb,double qeps)
  {
//--- save
   state.m_rstate.ia.Set(0,i);
   state.m_rstate.ia.Set(1,j);
   state.m_rstate.ia.Set(2,ns);
   state.m_rstate.ia.Set(3,info);
   state.m_rstate.ra.Set(0,c1);
   state.m_rstate.ra.Set(1,c2);
   state.m_rstate.ra.Set(2,intg);
   state.m_rstate.ra.Set(3,intk);
   state.m_rstate.ra.Set(4,inta);
   state.m_rstate.ra.Set(5,v);
   state.m_rstate.ra.Set(6,ta);
   state.m_rstate.ra.Set(7,tb);
   state.m_rstate.ra.Set(8,qeps);
  }
//+------------------------------------------------------------------+
//| Auxiliary function for AutoGKInternalIteration. Is a product to  |
//| get rid of the operator unconditional jump goto.                 |
//+------------------------------------------------------------------+
bool CAutoGK::Func_Internal_lbl_5(CAutoGKInternalState &state,
                                  int &i,int &j,int &ns,int &info,
                                  double &c1,double &c2,double &intg,
                                  double &intk,double &inta,double &v,
                                  double &ta,double &tb,double &qeps)
  {
//--- check
   if(i>state.m_n-1)
      return(Func_Internal_lbl_7(state,i,j,ns,info,c1,c2,intg,intk,inta,v,ta,tb,qeps));
//--- obtain F
   state.m_x=c1*state.m_qn[i]+c2;
   state.m_rstate.stage=0;
//--- Saving state
   Func_Internal_lbl_rcomm(state,i,j,ns,info,c1,c2,intg,intk,inta,v,ta,tb,qeps);
//--- return result
   return(true);
  }
//+------------------------------------------------------------------+
//| Auxiliary function for AutoGKInternalIteration. Is a product to  |
//| get rid of the operator unconditional jump goto.                 |
//+------------------------------------------------------------------+
bool CAutoGK::Func_Internal_lbl_7(CAutoGKInternalState &state,
                                  int &i,int &j,int &ns,int &info,
                                  double &c1,double &c2,double &intg,
                                  double &intk,double &inta,double &v,
                                  double &ta,double &tb,double &qeps)
  {
//--- calculation
   intk=intk*(state.m_b-state.m_a)*0.5;
   intg=intg*(state.m_b-state.m_a)*0.5;
   inta=inta*(state.m_b-state.m_a)*0.5;
   state.m_heap.Set(0,0,MathAbs(intg-intk));
   state.m_heap.Set(0,1,intk);
   state.m_heap.Set(0,2,inta);
   state.m_heap.Set(0,3,state.m_a);
   state.m_heap.Set(0,4,state.m_b);
   state.m_sumerr=state.m_heap[0][0];
   state.m_sumabs=MathAbs(inta);
//--- method iterations
   return(Func_Internal_lbl_14(state,i,j,ns,info,c1,c2,intg,intk,inta,v,ta,tb,qeps));
  }
//+------------------------------------------------------------------+
//| Auxiliary function for AutoGKInternalIteration. Is a product to  |
//| get rid of the operator unconditional jump goto.                 |
//+------------------------------------------------------------------+
bool CAutoGK::Func_Internal_lbl_8(CAutoGKInternalState &state,
                                  int &i,int &j,int &ns,int &info,
                                  double &c1,double &c2,double &intg,
                                  double &intk,double &inta,double &v,
                                  double &ta,double &tb,double &qeps)
  {
//--- check
   if(j>ns-1)
     {
      //--- method iterations
      return(Func_Internal_lbl_14(state,i,j,ns,info,c1,c2,intg,intk,inta,v,ta,tb,qeps));
     }
//--- calculation
   ta=state.m_a+j*(state.m_b-state.m_a)/ns;
   tb=state.m_a+(j+1)*(state.m_b-state.m_a)/ns;
   c1=0.5*(tb-ta);
   c2=0.5*(tb+ta);
   intg=0;
   intk=0;
   inta=0;
   i=0;
//--- function call, return result
   return(Func_Internal_lbl_11(state,i,j,ns,info,c1,c2,intg,intk,inta,v,ta,tb,qeps));
  }
//+------------------------------------------------------------------+
//| Auxiliary function for AutoGKInternalIteration. Is a product to  |
//| get rid of the operator unconditional jump goto.                 |
//+------------------------------------------------------------------+
bool CAutoGK::Func_Internal_lbl_11(CAutoGKInternalState &state,
                                   int &i,int &j,int &ns,int &info,
                                   double &c1,double &c2,double &intg,
                                   double &intk,double &inta,double &v,
                                   double &ta,double &tb,double &qeps)
  {
//--- check
   if(i>state.m_n-1)
     {
      //--- calculation
      intk=intk*(tb-ta)*0.5;
      intg=intg*(tb-ta)*0.5;
      inta=inta*(tb-ta)*0.5;
      state.m_heap.Set(j,0,MathAbs(intg-intk));
      state.m_heap.Set(j,1,intk);
      state.m_heap.Set(j,2,inta);
      state.m_heap.Set(j,3,ta);
      state.m_heap.Set(j,4,tb);
      state.m_sumerr=state.m_sumerr+state.m_heap[j][0];
      state.m_sumabs=state.m_sumabs+MathAbs(inta);
      j=j+1;
      //--- function call, return result
      return(Func_Internal_lbl_8(state,i,j,ns,info,c1,c2,intg,intk,inta,v,ta,tb,qeps));
     }
//--- obtain F
   state.m_x=c1*state.m_qn[i]+c2;
   state.m_rstate.stage=1;
//--- Saving state
   Func_Internal_lbl_rcomm(state,i,j,ns,info,c1,c2,intg,intk,inta,v,ta,tb,qeps);
//--- return result
   return(true);
  }
//+------------------------------------------------------------------+
//| Auxiliary function for AutoGKInternalIteration. Is a product to  |
//| get rid of the operator unconditional jump goto.                 |
//+------------------------------------------------------------------+
bool CAutoGK::Func_Internal_lbl_14(CAutoGKInternalState &state,
                                   int &i,int &j,int &ns,int &info,
                                   double &c1,double &c2,double &intg,
                                   double &intk,double &inta,double &v,
                                   double &ta,double &tb,double &qeps)
  {
//--- additional memory if needed
   if(state.m_heapused==state.m_heapsize)
      MHeapResize(state.m_heap,state.m_heapsize,4*state.m_heapsize,state.m_heapwidth);
//--- TODO: every 20 iterations recalculate errors/sums
   if(state.m_sumerr<=state.m_eps*state.m_sumabs || state.m_heapused>=m_maxsubintervals)
     {
      state.m_r=0;
      for(j=0; j<=state.m_heapused-1; j++)
         state.m_r=state.m_r+state.m_heap[j][1];
      //--- return result
      return(false);
     }
//--- Exclude interval with maximum absolute error
   MHeapPop(state.m_heap,state.m_heapused,state.m_heapwidth);
   state.m_sumerr=state.m_sumerr-state.m_heap[state.m_heapused-1][0];
   state.m_sumabs=state.m_sumabs-state.m_heap[state.m_heapused-1][2];
//--- Divide interval,create subintervals
   ta=state.m_heap[state.m_heapused-1][3];
   tb=state.m_heap[state.m_heapused-1][4];
   state.m_heap.Set(state.m_heapused-1,3,ta);
   state.m_heap.Set(state.m_heapused-1,4,0.5*(ta+tb));
   state.m_heap.Set(state.m_heapused,3,0.5*(ta+tb));
   state.m_heap.Set(state.m_heapused,4,tb);
   j=state.m_heapused-1;
//--- function call, return result
   return(Func_Internal_lbl_16(state,i,j,ns,info,c1,c2,intg,intk,inta,v,ta,tb,qeps));
  }
//+------------------------------------------------------------------+
//| Auxiliary function for AutoGKInternalIteration. Is a product to  |
//| get rid of the operator unconditional jump goto.                 |
//+------------------------------------------------------------------+
bool CAutoGK::Func_Internal_lbl_16(CAutoGKInternalState &state,
                                   int &i,int &j,int &ns,int &info,
                                   double &c1,double &c2,double &intg,
                                   double &intk,double &inta,double &v,
                                   double &ta,double &tb,double &qeps)
  {
//--- check
   if(j>state.m_heapused)
     {
      //--- function call
      MHeapPush(state.m_heap,state.m_heapused-1,state.m_heapwidth);
      //--- function call
      MHeapPush(state.m_heap,state.m_heapused,state.m_heapwidth);
      state.m_heapused=state.m_heapused+1;
      //--- method iterations
      return(Func_Internal_lbl_14(state,i,j,ns,info,c1,c2,intg,intk,inta,v,ta,tb,qeps));
     }
//--- calculation
   c1=0.5*(state.m_heap[j][4]-state.m_heap[j][3]);
   c2=0.5*(state.m_heap[j][4]+state.m_heap[j][3]);
   intg=0;
   intk=0;
   inta=0;
   i=0;
//--- function call, return result
   return(Func_Internal_lbl_19(state,i,j,ns,info,c1,c2,intg,intk,inta,v,ta,tb,qeps));
  }
//+------------------------------------------------------------------+
//| Auxiliary function for AutoGKInternalIteration. Is a product to  |
//| get rid of the operator unconditional jump goto.                 |
//+------------------------------------------------------------------+
bool CAutoGK::Func_Internal_lbl_19(CAutoGKInternalState &state,
                                   int &i,int &j,int &ns,int &info,
                                   double &c1,double &c2,double &intg,
                                   double &intk,double &inta,double &v,
                                   double &ta,double &tb,double &qeps)
  {
//--- check
   if(i>state.m_n-1)
     {
      //--- calculation
      intk=intk*(state.m_heap[j][4]-state.m_heap[j][3])*0.5;
      intg=intg*(state.m_heap[j][4]-state.m_heap[j][3])*0.5;
      inta=inta*(state.m_heap[j][4]-state.m_heap[j][3])*0.5;
      state.m_heap.Set(j,0,MathAbs(intg-intk));
      state.m_heap.Set(j,1,intk);
      state.m_heap.Set(j,2,inta);
      state.m_sumerr=state.m_sumerr+state.m_heap[j][0];
      state.m_sumabs=state.m_sumabs+state.m_heap[j][2];
      j=j+1;
      //--- function call, return result
      return(Func_Internal_lbl_16(state,i,j,ns,info,c1,c2,intg,intk,inta,v,ta,tb,qeps));
     }
//--- F(x)
   state.m_x=c1*state.m_qn[i]+c2;
   state.m_rstate.stage=2;
//--- Saving state
   Func_Internal_lbl_rcomm(state,i,j,ns,info,c1,c2,intg,intk,inta,v,ta,tb,qeps);
//--- return result
   return(true);
  }
//+------------------------------------------------------------------+
//| Iterative method                                                 |
//+------------------------------------------------------------------+
bool CAutoGK::AutoGKIteration(CAutoGKState &state)
  {
//--- create variables
   double s=0;
   double tmp=0;
   double eps=0;
   double a=0;
   double b=0;
   double x=0;
   double t=0;
   double alpha=0;
   double beta=0;
   double v1=0;
   double v2=0;
//--- This code initializes locals by:
//--- * random values determined during code
//---   generation - on first subroutine call
//--- * values from previous call - on subsequent calls
   if(state.m_rstate.stage>=0)
     {
      //--- initialization
      s=state.m_rstate.ra[0];
      tmp=state.m_rstate.ra[1];
      eps=state.m_rstate.ra[2];
      a=state.m_rstate.ra[3];
      b=state.m_rstate.ra[4];
      x=state.m_rstate.ra[5];
      t=state.m_rstate.ra[6];
      alpha=state.m_rstate.ra[7];
      beta=state.m_rstate.ra[8];
      v1=state.m_rstate.ra[9];
      v2=state.m_rstate.ra[10];
     }
   else
     {
      //--- initialization
      s=-983;
      tmp=-989;
      eps=-834;
      a=900;
      b=-287;
      x=364;
      t=214;
      alpha=-338;
      beta=-686;
      v1=912;
      v2=585;
     }
//--- check
   if(state.m_rstate.stage==0)
     {
      //--- change values
      state.m_needf=false;
      state.m_nfev=state.m_nfev+1;
      state.m_internalstate.m_f=state.m_f;
      //--- function call, return result
      return(Func_lbl_5(state,s,tmp,eps,a,b,x,t,alpha,beta,v1,v2));
     }
//--- check
   if(state.m_rstate.stage==1)
     {
      //--- change value
      state.m_needf=false;
      //--- check
      if(alpha!=0.0)
         state.m_internalstate.m_f=state.m_f*MathPow(x,-(alpha/(1+alpha)))/(1+alpha);
      else
         state.m_internalstate.m_f=state.m_f;
      state.m_nfev=state.m_nfev+1;
      //--- function call, return result
      return(Func_lbl_9(state,s,tmp,eps,a,b,x,t,alpha,beta,v1,v2));
     }
//--- check
   if(state.m_rstate.stage==2)
     {
      //--- change value
      state.m_needf=false;
      //--- check
      if(beta!=0.0)
         state.m_internalstate.m_f=state.m_f*MathPow(x,-(beta/(1+beta)))/(1+beta);
      else
         state.m_internalstate.m_f=state.m_f;
      state.m_nfev=state.m_nfev+1;
      //--- function call, return result
      return(Func_lbl_11(state,s,tmp,eps,a,b,x,t,alpha,beta,v1,v2));
     }
//--- Routine body
   eps=0;
   a=state.m_a;
   b=state.m_b;
   alpha=state.m_alpha;
   beta=state.m_beta;
   state.m_terminationtype=-1;
   state.m_nfev=0;
   state.m_nintervals=0;
//--- smooth function  at a finite interval
   if(state.m_wrappermode!=0)
     {
      //--- function with power-law singularities at the ends of a finite interval
      if(state.m_wrappermode!=1)
         return(false);
      //--- test coefficients
      if(alpha<=-1.0 || beta<=-1.0)
        {
         //--- change values
         state.m_terminationtype=-1;
         state.m_v=0;
         //--- return result
         return(false);
        }
      //--- special cases
      if(a==b)
        {
         //--- change values
         state.m_terminationtype=1;
         state.m_v=0;
         //--- return result
         return(false);
        }
      //--- reduction to general form
      if(a<b)
         s=1;
      else
        {
         //--- change values
         s=-1;
         tmp=a;
         a=b;
         b=tmp;
         tmp=alpha;
         alpha=beta;
         beta=tmp;
        }
      //--- change values
      alpha=MathMin(alpha,0);
      beta=MathMin(beta,0);
      //--- first,integrate left half of [a,b]:
      //---     integral(f(x)dx,a,(b+a)/2)=
      //---     =1/(1+alpha) * integral(t^(-alpha/(1+alpha))*f(a+t^(1/(1+alpha)))dt,0,(0.5*(b-a))^(1+alpha))
      AutoGKInternalPrepare(0,MathPow(0.5*(b-a),1+alpha),eps,state.m_xwidth,state.m_internalstate);
      //--- function call, return result
      return(Func_lbl_9(state,s,tmp,eps,a,b,x,t,alpha,beta,v1,v2));
     }
//--- special case
   if(a==b)
     {
      //--- change values
      state.m_terminationtype=1;
      state.m_v=0;
      //--- return result
      return(false);
     }
//--- general case
   AutoGKInternalPrepare(a,b,eps,state.m_xwidth,state.m_internalstate);
//--- function call, return result
   return(Func_lbl_5(state,s,tmp,eps,a,b,x,t,alpha,beta,v1,v2));
  }
//+------------------------------------------------------------------+
//| Auxiliary function for AutoGKIteration. Is a product to get rid  |
//| of the operator unconditional jump goto.                         |
//+------------------------------------------------------------------+
void CAutoGK::Func_lbl_rcomm(CAutoGKState &state,double s,double tmp,
                             double eps,double a,double b,double x,
                             double t,double alpha,double beta,
                             double v1,double v2)
  {
//--- save
   state.m_rstate.ra.Set(0,s);
   state.m_rstate.ra.Set(1,tmp);
   state.m_rstate.ra.Set(2,eps);
   state.m_rstate.ra.Set(3,a);
   state.m_rstate.ra.Set(4,b);
   state.m_rstate.ra.Set(5,x);
   state.m_rstate.ra.Set(6,t);
   state.m_rstate.ra.Set(7,alpha);
   state.m_rstate.ra.Set(8,beta);
   state.m_rstate.ra.Set(9,v1);
   state.m_rstate.ra.Set(10,v2);
  }
//+------------------------------------------------------------------+
//| Auxiliary function for AutoGKIteration. Is a product to get rid  |
//| of the operator unconditional jump goto.                         |
//+------------------------------------------------------------------+
bool CAutoGK::Func_lbl_5(CAutoGKState &state,double &s,double &tmp,
                         double &eps,double &a,double &b,double &x,
                         double &t,double &alpha,double &beta,
                         double &v1,double &v2)
  {
//--- check
   if(!AutoGKInternalIteration(state.m_internalstate))
     {
      //--- change values
      state.m_v=state.m_internalstate.m_r;
      state.m_terminationtype=state.m_internalstate.m_info;
      state.m_nintervals=state.m_internalstate.m_heapused;
      //--- return result
      return(false);
     }
//--- change values
   x=state.m_internalstate.m_x;
   state.m_x=x;
   state.m_xminusa=x-a;
   state.m_bminusx=b-x;
   state.m_needf=true;
   state.m_rstate.stage=0;
//--- Saving state
   Func_lbl_rcomm(state,s,tmp,eps,a,b,x,t,alpha,beta,v1,v2);
//--- return result
   return(true);
  }
//+------------------------------------------------------------------+
//| Auxiliary function for AutoGKIteration. Is a product to get rid  |
//| of the operator unconditional jump goto.                         |
//+------------------------------------------------------------------+
bool CAutoGK::Func_lbl_9(CAutoGKState &state,double &s,double &tmp,
                         double &eps,double &a,double &b,double &x,
                         double &t,double &alpha,double &beta,
                         double &v1,double &v2)
  {
//--- check
   if(!AutoGKInternalIteration(state.m_internalstate))
     {
      v1=state.m_internalstate.m_r;
      state.m_nintervals=state.m_nintervals+state.m_internalstate.m_heapused;
      //--- then,integrate right half of [a,b]:
      //---     integral(f(x)dx,(b+a)/2,b)=
      //---     =1/(1+beta) * integral(t^(-beta/(1+beta))*f(b-t^(1/(1+beta)))dt,0,(0.5*(b-a))^(1+beta))
      AutoGKInternalPrepare(0,MathPow(0.5*(b-a),1+beta),eps,state.m_xwidth,state.m_internalstate);
      //--- function call, return result
      return(Func_lbl_11(state,s,tmp,eps,a,b,x,t,alpha,beta,v1,v2));
     }
//--- Fill State.X,State.XMinusA,State.BMinusX.
//--- Latter two are filled correctly even if B<A.
   x=state.m_internalstate.m_x;
   t=MathPow(x,1/(1+alpha));
   state.m_x=a+t;
//--- check
   if(s>0.0)
     {
      state.m_xminusa=t;
      state.m_bminusx=b-(a+t);
     }
   else
     {
      state.m_xminusa=a+t-b;
      state.m_bminusx=-t;
     }
//--- change values
   state.m_needf=true;
   state.m_rstate.stage=1;
//--- Saving state
   Func_lbl_rcomm(state,s,tmp,eps,a,b,x,t,alpha,beta,v1,v2);
//--- return result
   return(true);
  }
//+------------------------------------------------------------------+
//| Auxiliary function for AutoGKIteration. Is a product to get rid  |
//| of the operator unconditional jump goto.                         |
//+------------------------------------------------------------------+
bool CAutoGK::Func_lbl_11(CAutoGKState &state,double &s,double &tmp,
                          double &eps,double &a,double &b,double &x,
                          double &t,double &alpha,double &beta,
                          double &v1,double &v2)
  {
//--- check
   if(!AutoGKInternalIteration(state.m_internalstate))
     {
      //--- change values
      v2=state.m_internalstate.m_r;
      state.m_nintervals=state.m_nintervals+state.m_internalstate.m_heapused;
      //--- final result
      state.m_v=s*(v1+v2);
      state.m_terminationtype=1;
      //--- return result
      return(false);
     }
//--- Fill State.X,State.XMinusA,State.BMinusX.
//--- Latter two are filled correctly (X-A,B-X) even if B<A.
   x=state.m_internalstate.m_x;
   t=MathPow(x,1/(1+beta));
   state.m_x=b-t;
//--- check
   if(s>0.0)
     {
      state.m_xminusa=b-t-a;
      state.m_bminusx=t;
     }
   else
     {
      state.m_xminusa=-t;
      state.m_bminusx=a-(b-t);
     }
//--- change values
   state.m_needf=true;
   state.m_rstate.stage=2;
//--- Saving state
   Func_lbl_rcomm(state,s,tmp,eps,a,b,x,t,alpha,beta,v1,v2);
//--- return result
   return(true);
  }
//+------------------------------------------------------------------+
